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Numerical  Model  for  Forecasting  Ice  Conditions 
on  the  Ohio  River 


HUNG  TAG  SHEN,  GORANKA  BJEDOV,  STEVEN  F.  DALY  AND  A.  M.  WASANTHA  LAL 


INTRODUCTION 

During  the  winter  in  areas  of  higher  latitude,  ice  can  form  in  rivers.  The  presence  of  ice  can 
significantly  influence  flow  and  navigation  conditions  in  the  river.  The  ability  to  forecast  river  ice 
conditions  is  therefore  of  great  importance  in  planning  winter  flow  regulation  and  navigation 
operations. 

The  ice  condition  in  a  river  not  only  influences,  but  also  interacts  with,  its  flow  condition.  Numerous 
computer  models  simulate  unsteady  flow  in  rivers  under  open-water  conditions  (Mah.nood  and 
Yevjevich  1975,  Cunge  et  al.  1980,  Fread  1985).  Yapa  and  Shen  (1986)  developed  an  unsteady  flow 
model  for  ice-covered  rivers  by  including  the  effects  of  hydraulic  resistance  of  the  ice  cover.  In  the 
model  the  transport  of  moving  ice  in  the  river  and  the  ice  cover  formation  were  not  considered, 
although  the  thermal  growth  and  decay  of  the  ice  cover  were  simulated  using  a  degree-day  model.  A 
few  models  are  capable  of  simulating  river  ice  conditions  (Marcotte  1981,  Michel  and  Drouin  1981, 
Petryk  1981,  Calkins  1984,  Shen  and  Yapa  1984).ExceptforthatofShenand  Yapa(1984),all  of  these 
models  use  a  backwater  computation  and  ignore  the  distribution  and  transport  of  ice  along  the  river. 
For  long  rivers  subject  to  repeated  freezing  and  melting  during  a  winter,  a  river  ice  model  must  simulate 
the  transport  process  correctly.  In  this  study  a  numerical  model  called  RICEOH  for  simulating  river 
ice  and  flow  conditions  is  developed  for  the  Ohio  River  system  between  Pittsburgh,  Pennsylvania,  and 
Meldahl,  Ohio.  With  weather  forecasts  as  input,  the  model  can  be  used  to  forecast  river  ice  conditions. 

The  model  is  based  on  a  single-channel  river  ice  model  recently  developed  by  Lai  (1988).  In  the 
hydraulics  computation  the  one-dimensional  unsteady  flow  model  developed  by  Chen  and  Simons 
(1975)  and  implemented  for  the  Ohio  River  under  open-water  conditions  (Johnson  1982)  is  used. 
Modifications  are  made  in  the  unsteady  flow  model  to  include  the  hydraulic  resistance  of  the  ice  cover. 
In  the  ice  model  the  water  temperature  and  ice  concentration  distributions  along  the  river  are  calculated 
by  a  Lagrangian-Eulerian  scheme.  The  formation  of  the  ice  cover  is  modeled  using  existing 
equilibrium  ice  jam  theories  (Pariset  and  Hausser  1961).  The  thermal  growth  and  decay  of  the  ice  cover 
are  computed  based  on  quasi-steady  thermal  conduction  in  the  ice  cover  considering  heat  exchanges 
between  the  atmosphere,  the  ice  cover  and  the  underlying  river  water  (Shen  and  Lai  1986). 
Modifications  are  made  on  the  single-channel  model  of  Lai  (1988)  to  make  the  model  applicable  to 
a  river  system  with  dendritic  tributaries. 


HYDRAULIC  ANALYSIS 


The  continuity  and  momentum  equations  for  a  river  with  a  floating  ice  cover,  as  shown  in  Figure 
1 ,  are  given  as 


dQ  ^  ^ 
dx  dt 


qi  =  0 


P 


dt 


+  p 


A  dx 


^  ^  +  (Pi'^i  +  Pb'tb)  -P<7ivi  =  0 

dx]  dx 


where  Q  = 
A  = 

X  = 
t  = 

Qi  = 
H  = 

y  = 

g  = 
= 

p  = 


discharge 
flow  area 

horizontal  distance  along  the  channel 
time 

lateral  inflow  to  the  channel 
water  level  {H  =  +  d^+ 

depth  of  the  water  (y-d^+ 
acceleration  due  to  gravity 
bed  elevation 
density  of  water 
depth  of  the  flow 

wetted  perimeter  formed  by  the  channel  bed 
wetted  perimeter  formed  by  the  ice  cover 
shear  stress  at  the  channel  bottom 
shear  stress  at  the  ice/water  interface 

lateral  inflow  velocity  component  in  the  main  stream  direction 
submerged  thickness  of  the  ice  cover. 


(1) 

(2) 


Figure  1.  Channel  flow  with  a  floating  ice  cover. 


Since  eq  1  and  2  do  not  in  general  possess  analytical  solutions,  one  must  rely  on  numerical 
techniques  to  solve  them.  In  the  present  study  the  computer  model  developed  for  the  Ohio  River  by 
Chen  and  Simons  (1975)  and  Johnson  (1982)  is  adopted  to  solve  eq  1  and  2,  with  modifications  for 
the  terms  related  to  the  effect  of  the  ice  cover.  The  numerical  method  used  is  a  linear  implicit  finite- 
difference  method  with  a  double-sweep  algorithm. 
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Figure  2.  Four-point  implicit  scheme. 


Finite-difTerence  formulation 

Figure  2  and  eq  3  show  the  discretization  of  dependent  variables  and  its  derivatives  in  time  and  space: 


=  1 
4 
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i 

1  r-n+i 

r/l+l  .  rlt  rttX 
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/ 
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where  At  is  the  time  step  used  in  computations  and  Ax  is  the  distance  between  two  neighboring  points. 
We  can  rearrange  the  one-dimensional  flow  equations  in  the  following  form: 


dt  dx 


and 


(4) 


2u^  -  u-B^  +gA^  =gA  (5o-5f)+  ^,V|  + 

av  av  '  "  ^  lav  1^=,, 


(5) 


where  6  is  the  channel  width,  u  =  QIA  is  the  averaged  cross-sectional  velocity  and  5^  and  5|- represent 
the  channel-bed  slope  and  energy  slope,  respectively.  According  toeq  3  these  equations  can  be  written 
in  finite-difference  form  as 


1 

2A.V 


[(er. ,  -  Q!)  *  (or;l  -  or*')]  *  b".,  ^  |[k)"*'  -  k)*] 
-Kir.,]) =0 
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-(«-b)4  ^  [K)r. .  -K)!'  +  (^4v):':'i -{d.yr\ 


A  1  r  /I  /I  .  /I + 1  /I + 1 

n  + ' 


(7) 


Ineq6and7,/,+i  =  ]^(f\  +/"+])  and/17  =  ( 3/1 /a.v )  ^  If  the  spatial  grid  has  different 

lengths,  then  Ax  in  the  above  equations  has  to  be  replaced  by  AVj. 

The  presence  of  an  ice  cover  can  affect  the  frictional  slope.  For  a  reach  as  shown  in  Figure  3,  the 
friction  slope  is  expressed  as 


“  “o.i  ('^f)o,i  “i,i 


(8) 


where  a  j  and  a  j  ^  are  the  length  fractions  of  open  water  and  ice  cover,  respectively,  in  the  length 
element  A/j  corresponding  to  node  i. 

The  friction  slope  of  the  open-water  portion  can  be  calculated  as 


i  ~ 


(^o,  i  Qo,  i) 


(l.486Ao,i/?oj3f 


(9) 


where  the  hydraulic  radius  is  defined  as 
^  o.  i 


^  o,  i  “ 


PoA 


(10) 


The  subscript  o  represents  the  open-water  conditions  with  no  ice  cover  effects,  and 
Pq  i  "  wetted  perimeter  of  node  i 
Ag  j  =  flow  area  at  node  / 

/tg  j  =  Manning’s  roughness  coefficient  of  the  channel  bed  for  node  i 
=  velocity  at  node  /. 

The  friction  slope  of  an  ice-covered  reach  is  defined  as 
(”c,  i  Qi.  i)' 


H.,  = 


(11) 


i-  1 


i  "  I 


Al: 


Ax. 


Figures.  Definition  sketch  for  the 
'  *  ’  friction  slope  calculation. 
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where  the  composite  roughness  coefficient  for  node  /  is 


2/3 


(12) 


where  j  is  the  undersurface  roughness  of  the  ice  cover  for  node  /. 

The  length  fraction  a  j  j  is  calculated  as 


a.  .  =  2 


ttiAa:;  +  a,-  1  AAri_  \ 


Ajtj  +  AjCi_  I 

where  coefficients  a  are  defined  in  Figure  3: 


(13) 


=X(i-l) 

for 

X(i-1)<0.5 

=  0.5 

for 

X(i-1)>0.5 

«i  =0 

for 

X(i)  <  0.5 

=  X(i>-0.5 

for 

X(i)>0.5 

and  X;  and  X- 


i-i 


AXj  and  Arj  j 
i ,  1-1  and  i+l 


=  fractions  of  the  length  covered  by  ice  in  reaches  i  and  i-l ,  namely  the  length 
of  the  ice  cover  divided  by  the  length  of  the  corresponding  reach 
=  submerged  thicknesses  of  ice  cover  in  reaches  i  and  i-l ,  respectively 
=  lengths  of  reaches  i  and  i-l,  respectively 
=  nodal  point  indices. 


The  length  fraction  Ax^  j  is 
Ajio.i  =  1  -  Axj  j. 

The  velocity  for  the  fully  ice-covered  condition  is  expressed  as 


where  the  flow  area  Aj  j  is  given  as 

i,  i  ~  i  —  Tj  /isu,  i  • 


(14) 


(15) 


(16) 


The  equivalent  submerged  ice  thickness  for  node  i  is  defined  as 

'  Ot  i  _  1  Ax  i  -  1  su,  i  -  1  ®  i  ^  ^  i  ^  su.  i 

= - 

■  ttj.  1  Aj:i_  I  +  ajAJCj 

and  the  hydraulic  radius  for  the  ice-covered  case  is 


Ri.i  = 


Pi,i 


(17) 


(18) 


where  the  wetted  perimeter  in  the  ice-covered  case  =  i>  ihc  top  width 

of  the  channel  at  node  /.  The  equivalent  ice  cover  roughness  coefficient  at  node  i  is  calculated  from 

ai_|Ajri_,/-j_i  -t-ttiAxir; 

Wj  i  = -  (IV) 

ai_  I  Axj_  1  +  OjAxi 
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where  r-^  and  /•j_,  are  ice  cover  roughness  coefficients  for  reaches  i  and  /-I ,  respectively. 

For  the  boundary  node  of  each  tributary,  eq  17-19  are  modified  accordingly.  For  the  upstream 
boundary  node  a  ^  j  =  2(X^  -  0.5)  forX,  >0.5  and  0  otherwise,  and  the  equivalent  submerged  thickness 
and  roughness  of  ice  cover  are  A  j  =  A  su ,  i- 1  i  ~  downstream  boundary  node  these 

parameters  are  defined  as  a  •  j  =  2Xj_|  forX|_j  <  0.5  and  1  otherwise,  and  h  |  =  ^  su,i-i  ^nd  Wj  [  = 
ri_,,  where  the  subscript  /  denotes  the  last  downstream  node  on  the  branch. 

To  assure  the  stability  of  the  numerical  scheme  (Strelkoff  1970),  the  friction  slo|)e  Sf  is  taken  on 
t  time  level  with  a  first-order  Taylor  expansion.  A  suitable  expression  for  S  f  is 


Q'l) 


,  <  n  +  l 

whereoCj  jisdefinedineq  13.Sincethelengthfractionofanicecovera.  .  is  unknown  at  the  present 
time  step,  it  is  necessary  to  approximate  the  last  term  in  eq  20.  The  approximation  that  was  used  in  this 
model  is 


-n  +  l 


-«-l 


a.  .  -a  =a.  .-a.. 
1.1  1,1  1,1  1,1 


Using  eq  8-19,  terms  in  eq  20  can  be  expressed  as  eq  21-23  considering  n  as  a  function  of  Q  and  y: 


3Sf  _ 


95  f  _ 

9y 


2(5 f)o  (J-  +  -L-  a;  +  2(Sf)i  ( L  +  X  a.' 

Iq  «o  90/  \Q  «c  dQl  ' 
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_l( 

Ao  '  3  3  dy  I  ”o  9y  . 


(21) 


and 


1 

"c  9>’. 


(22) 


^  =(5f)i-(5f)o.  (23) 

Substituting  eq  10-1 3  into  eq  6  and  7  yields 

-  00i  +  (Ci)i  yj  +  00i+i  +  (C i)i  yi+  1  =(C2)i  (24) 

and 
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where  0  =  ^ 
Ajt 


(c,),  = 


2  2  J 
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All  coefficients  in  eq  24  and  25  are  evaluated  at  the  time  level  /”  and  therefore  are  known. 

There  are  four  unknowns  in  eq  24  and  25  at  the  time  level  r However,  two  unknowns  are 
common  for  any  two  neighboring  grid  points.  Consequently  the  (/-I)  pairs  of  equations  contain  2/ 
unknowns.  Therefore,  we  need  two  additional  equations  that  can  be  established  from  boundary 
conditions  to  complete  this  system. 

At  the  upstream  boundary  it  is  assumed  that  the  flood  hydrograph  supplies  the  flow  depth  or 
discharge  as  a  function  of  time.  This  relation  can  be  written  as 

C9  (2i  +  Cjo  yi  =  Cii .  (26) 

At  the  downstream  boundary  the  rating  curve  is  expressed  in  segmentized  form  as 

C12  Gi  +  C|3  yi  =  Ci4 .  (27) 

Equations  24-27  form  a  system  of  2/  linear  algebraic  equations  with  2/  unknowns.  Any  standard 
method  can  be  used  for  its  solution.  In  this  model  the  double-sweep  algorithm  was  used. 

Double-sweep  algorithm 

If  we  rearrange  eq  24  and  25  by  subtracting  from  them  the  same  set  of  equations  from  the  previous 


time  step,  we  get 

A/y  Ay  i+i  +  Bj  A(3i^.i  =  C  j  Avj  +  DiAQ\  +  G\  (28) 

W^Ayi^,  +  =  c'.Ay,+DlAQi+G'.  (29) 

where  Ay  i  ^y"^*  -  y”  and  AQ^  -  -Oi”- 

Assume  now  that  there  is  a  linear  relationship  of  the  type 

A(2i  -  EiAyi  +Fi  (30) 

for  point  /.  It  can  be  easily  proven  that  an  analogous  linear  relation.ship  exists  for  the  next  point  (+ 1  and 
a  recurrence  relationship,  eq  3 1 ,  can  be  obtained  (Liggett  and  Cunge  1 975): 

'^Qi  +  l  =  B^i  +  l'^>’i  +  l  +  A^i  +  l  •  (31) 

Coefficients  and  can  be  computed  if  Fj  and  F-  are  known.  If  we  now  substitute  eq  30  in  eq 
28,  we  can  get 

Ayi  =  FiAyj+,  -HMiAGi  +  l  -t-Ai  (32) 


where  L  -  - CIA - 

Ci-i-  D.Fj 
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Afj  =  _ £j _ 

C  i  +  DyE\ 

yy  _  G\  +  D\F\ 

'  Ci+Di£i' 

Equation  31  permits  the  computation  of  A^j  when  the  increments  Ajj^j  and  AO^^I  are  known. 
Therefore,  it  is  possible  to  compute  and  for  all  points  of  a  given  reach.  This  method  is  fully 
explained  in  the  block  diagram  of  Figure  4.  The  interaction  between  the  main  river  and  a  tributary  and 
the  treatment  of  locks  and  dams  are  explained  in  detail  by  Johnson  (1982).  The  practical  meaning  of 
the  double-sweep  algorithm  is  that  the  number  of  elementary  operations  needed  to  solve  the  system 
is  proportional  to  the  number  of  points  I.  Standard  methods,  such  as  matrix  inversion,  have  the  number 
of  operations  proportional  to  (Burden  and  Faires  1986). 


i  =  1 


Compute  E  j ,  Fj  From 
Boundary  Condition  at 
Point  i  =  1 


Compute  Coemcients  L| ,  Mj  and  N| 
Compute  Coefficients  E  j  +  ,  and  Fj 
Store  Lj,  M  j.  N  |,  Ei  +  1  and  F  , 


Compute  Dz  |  from  Boundary 
Condition  at  Point  i  +  1  =  I 


Compute  DQ| ,  Z|  and  Q| 


Compute  Dzjand  OOj  and 
Qi  and  z  j  at  New  Time  Step 


Figure  4.  Block  diagram  for  the  double-sweep  algorithm. 
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After  the  flow  conditions  at  the  time  level  are  computed  for  each  node  by  solving  eq  4  and  5, 
the  temperature  and  ice  conditions  are  calculated  along  the  river.  The  procedures  for  thermal  and  ice 
computations  are  described  in  the  following  section. 


SIMULATION  OF  THERMAL 
AND  ICE  CONDITIONS 

The  present  model  adopts  the  model  developed  by  Lai  ( 1 988)  to  simulate  thermal  and  ice  conditions, 
with  some  modifications. 

Water  temperature 

and  ice  discharge  distributions 

In  a  well-mixed  river  the  governing  equation  for  the  distribution  of  water  temperature  along  the 

river  is 

A(pCpAr^)  +  |-(epCprJ  =  A(A£,pCp^)  +  B<i>  +  i7,pCp(r,-TJ  (33) 
at  ax  dv  \  dx  j 

where  A  =  flow  area 
B  =  river  width 
p  =  density  of  water 
Cp  =  specific  heat  of  water 

=  longitudinal  dispersion  coefficient 
<t>  =  heat  flux  per  unit  surface  area  of  the  river 
r,  =  temperature  of  the  lateral  inflow. 

On  the  river  surface  the  energy  flux  consists  of  solar  or  short-wave  radiation,  long- wave  radiation, 
heat  transfer  due  to  evaporation  or  condensation,  sensible  heat  transfer  due  to  conduction,  and  heat 
transfer  due  to  precipitation.  A  simplified  linear  model  is  used  in  this  study.  The  net  heat  exchange  at 
the  free  surface  is  expressed  as  (Ashton  1986) 


^  —  ^*wa(^a  ~  ^w) 


(34) 


where  =  water  temperature 
=  air  temperature 

/iwa  =  heat  exchange  coefficient  at  the  water/air  interface. 

Linear  models  cannot  accurately  describe  the  heat  exchange  process.  However,  for  rivers  where 
extensive  weather  data  are  not  available,  they  provide  sufficiently  accurate  prediction  of  water 
temperature. 

When  a  river  is  covered  with  ice,  the  turbulent  heat  exchange  between  the  river  water  and  the  ice 
cover  is  described  by  (Ashton  1986) 

Owi  =  ^Iwi  (fw  ~  fm)  (35) 


where  is  the  melting  temperature  (0°C)  and  is  the  turbulent  heat  exchange  coefficient  given  as 


h  -  r  U- 

'•  W1  “  W1 


0.8 


0.2 


(36) 
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where  the  coefficient  C^^,j  approximately  equals  1622  W  m~^-^  °C“'.  Laboratoiy  and  field 
investigations  indicate  that  may  increase  with  the  resistance  of  the  ice  cover. 

When  the  water  temperature  drops  to  the  freezing  point  T^,  frazil  ice  production  starts  in  a  river.  The 
equation  for  frazil  ice  transport  can  be  written  as 

l-(piLiCi)  +  A(2  piLiCi)  =  +  Bh^.ATa-Tf)  (37) 

3/  dx  ^  '  dx\  dx  I  ' 

where  Cj  =  ice  concentration 

L-  =  latent  heat  of  fusion 
pj  =  density  of  ice. 

Since  eq  37  is  in  the  same  form  as  eq  33,  the  soluticoi  of  eq  33  when  <  0°C  can  be  used  to  determine 

the  ice  concentration  by  setting  =  Tf  in  eq  34  and  letting  Cj  =  pCpf  w  /  Pi^  i  in  the  solution. 

If  we  neglect  the  dispersion  term,  eq  33  can  be  reduced  to  eq  38  by  using  the  continuity  equation 

(eq  1): 

|-  (pCpT-w)  +  Q^  (pCpT-w)  =  fiO  +  ^1  P  Cp(Ti  -  7^)  .  (38) 

In  Lagrangian  form  eq  38  becomes 

DT^  ^  O  ^ 

Dt  pCpd^  A 

In  the  present  model  eq  39  is  solved  using  a  Lagrangian-Eulerian  scheme.  In  this  scheme  parcels 
of  water  with  known  water  temperature  or  ice  concentration  at  time  f  are  followed  along  the  river  to 
obtain  the  temperature  or  ice  concentration  distribution  at  t”+At.  As  shown  in  Figure  5,  a  water  parcel 
located  at  Xj  at  t"  will  move  to  a  new  location  at  time  t"+At.  The  x  coordinate  of  this  new  location  is 
given  by 

7-1  /  7-1  \ 

Si=Xi+  51  Z  S'k 

t:=l  \  k  =l  } 

where  Sj  =  new  position  of  the  panicle 

5t|^  =  “  travel  time  of  the  parcel  in  reach 

y-1  =  last  node  passed  by  the  moving  parcel. 


Figure  5.  The  L(igran_^ian-Eulenan  scheme. 
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The  water  temperature  of  the  parcel  can  be  obtained  by  numerically  integrating  eq  39; 


^  =  (rw)!'  -  I  5,^ 

/'  =  !  \P^n«u/k  ^ 


where  '  is  the  temperature  of  the  water  parcel,  which  was  originally  located  at  .Vj,  and  is 

the  mean  depth  of  reach  k.  At  the  end  of  each  time  step,  the  water  temperature  at  every  grid  point  is 
obtained  by  linearly  interpolating  from  the  closest  Lagrangian  parcels  upstream  and  downstream  of 
the  grid  point.  These  interpolated  values  will  be  used  as  starting  points  for  the  next  time  step.  An  outline 
of  the  procedure  is  given  in  Figure  5. 

When  A/  is  large,  the  water  parcel  originated  from  the  upstream  boundary  at  time  t "  can  travel  over 
several  grid  points.  There  will  be  no  Lagrangian  parcels  located  in  the  upstream  portion  of  the  river 
for  interpolation.  Water  temperatures  of  nodal  points  in  this  region  are  determined  by  tracing  back  in 
time  to  the  upstream  boundary,  where  the  original  temperatures  of  water  parcels  are  known  at  all  times. 

When  calculating  the  water  temperature  for  the  first  node  in  the  main  stream  downstream  of  a 
junction,  instantaneous  mixing  of  water  from  the  tributary  with  that  of  the  main  stream  is  assumed; 


<I> 


pC  p^fwj 


A 


At  -  X 

k  =  \ 


(41) 


T 


,  _  T wl(2wl  +  T vjlQwl 
w3 - - 

Cwl  +  Qw2 


(42) 


where  Q^^  and  Q^2 

T^^  andr^2 


=  discharges  in  the  main  stream  and  the  tributary,  respectively 
=  water  temperatures  in  the  main  stream  and  the  tributary  immediately 
upstream  of  the  junction 

=  water  temperature  in  the  main  stream  immediately  downstream 
of  the  junction. 


Ice  cover  formation 

The  initial  ice  cover  formation  is  a  process  governed  by  channel  geometry,  hydraulic  conditions, 
surface  ice  supply  and  ice  properties.  Depending  on  these  conditions,  different  modes  of  ice  cover 
formation  can  occur  on  the  river.  The  present  model  is  capable  of  simulating  ice  cover  formation  by 
particle  juxtaposition,  hydraulic  thickening  (narrow  jam  formation)  and  mechanical  thickening  (wide 
jam  formation).  After  the  ice  cover  is  initiated,  it  will  progress  upstream,  depending  on  the  thickness 
of  the  ice  cover  and  the  supply  of  ice  suspended  in  water.  When  the  ice  cover  progresses  to  a  junction, 
it  can  progress  into  both  the  tributary  and  the  main  stream.  During  this  process  the  ice  cover  can  change 
its  thickness  due  to  thermal  growth  or  decay  and  deposition  or  erosion  of  frazil  ice  underneath  the  ice 
cover. 

In  regions  with  low  flow  velocity,  ice  covers  can  form  by  simple  juxtaposition  of  ice  floes.  The 
stability  condition  for  incoming  ice  floes  at  the  leading  edge  is  given  by  (Ashton  1974) 
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where  =  critical  velocity  upstream  of  the  leading  edge  for  undertuming  and  submergence 
T;  =  thickness  of  the  ice  floe 
y  =  upstream  flow  depth  at  the  leading  edge. 

Another  expression  has  been  obtained  earlier  by  Pariset  and  Hausser  (1961): 


F 


rc 


Ui/\ 


where 


=  form  factor,  which  varies  between  0.6  and  1 .3 
=  length  of  the  ice  floe 
=  porosity  of  the  ice  floe. 


(44) 


More  recently,  Daly  and  Axelson  (1990)  presented  a  refined  formulation  for  the  juxtaposition 
phenomenon.  If  the  criterion  for  juxtaposition  is  satisfied,  the  ice  cover  will  progress  upstream  with 
the  same  thickness  as  the  incoming  ice  floes. 

Field  observations  (Kivisild  1959)  indicated  that  F^^  can  vary  from  0.05  to  0. 10,  depending  on  the 
floe  characteristics.  Since  analytical  formulas  forF^^  require  the  geometry  of  the  ice  floe,  it  is  difficult 
to  apply  them  to  field  problems.  In  the  present  model  the  value  of  the  critical  Froude  number  F^^  is 
considered  to  be  an  input  parameter  that  specifies  the  limiting  condition  for  the  juxtaposition  mode. 
Floe  thickness  and  porosity  e,  which  are  needed  in  determining  the  rate  of  progression,  are  also 
specified. 

When  the  Froude  number  becomes  greater  than  the  critical  Froude  number,  the  ice  cover  will 
progress  as  a  narrow  jam  or  by  hydraulic  thickening.  Incoming  floes  will  submerge  and  become 
deposited  on  the  underside  of  the  ice  cover  to  a  thickness  governed  by  the  flow  condition. 

To  obtain  the  thickness  of  an  ice  cover  in  the  hydraulic  thickening  mode,  it  is  necessary  to  consider 
the  interaction  between  the  icecover  formation  andthe  flow  condition.  Ifwe  apply  the  energy  equation 
between  sections  1  and  3  on  Figure  6,  we  get 


.V| 


1 


-  ^w3  +  —  + 

2  a 


El 

pg  ■ 


(45) 


Fif’iire  6.  Definition  sketch  for  the  narrow  jam  formulation. 
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Figure  7.  Definition  sketch  for  an  ice-covered  reach. 


A  no-spill  condition  at  section  2  gives 
2 

y  1  +  ^  <  dw3  +  h  ■  (46) 

Using  a  hydrostatic  condition  the  pressure  underneath  the  ice  cover  at  section  3,  p^,  can  be  obtained 
from 

P3  =  pg/i  -  (1 -e)(p- Pi)gA  .  (47) 

By  substituting  eq  45  and  46  into  eq  47  we  obtain  the  first  equation  for  A  and  h: 

F,(/l,h)=Q2_2g^2ji_^j;,=0.  (48) 

Since  the  water  level  interacts  with  the  thickening  of  the  ice  cover,  an  additional  equation  is  needed 
for  i4  and  h.  An  equation  for  hydraulic  conditions  is  obtained  by  applying  the  energy  equation  between 
points  1  and  2  of  Figure  7: 

2  2 

H,  +-2—  =  //2  +  -^+ SfA;r  .  (49) 

2gA f  IgA  i 

If  we  express  the  net  flow  area  as 

A  =A'-  B  ^h  (50) 

P 

where  A'  is  the  known  flow  area  corresponding  to  the  depth  of  the  flow  y,  and  substitute  for  Sf  in  the 
eq  49,  we  can  obtain  a  hydraulics  backwater  equation: 

F2(-4,*)  =  0.  (51) 

Equations  48  and  51  are  solved  simultaneously  for  A  and  h  using  a  Newton-Raphson  procedure. 
In  the  narrow  jam  mode,  there  is  a  limiting  Froude  number  F,^^  beyond  which  ice  cover  caiuiot 
progress.  The  value  of  which  typically  equals  0.09,  is  considered  to  be  an  input  to  this  model. 

If  the  net  streamwise  force  exceeds  the  internal  resistance  of  the  ice  cover  formed  by  hydraulic 
thickening,  the  cover  will  collapse  and  thicken  until  an  equilibrium  thickness  is  reached.  This  process 
is  commonly  known  as  mechanical  thickening  of  an  ice  cover,  or  “shoving,"  and  accumulations 
formed  in  this  manner  are  often  called  wide  river  jams.  When  shoving  occurs  on  the  river,  a  relatively 
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long  reach  of  ice  cover  will  collapse,  and  the  leading  edge  will  move  downstream.  Based  on  the 
analyses  of  Pariset  and  Hausser  (1961)  and  Uzuner  and  Kennedy  (1976),  an  equation  for  the 
equilibrium  thickness  of  wide  jam  t  for  a  steady  uniform  flow  can  be  obtained  as  follows: 

where /j  and /j,  are  Darcy-Weisbach  friction  factors  related  to  the  ice  cover  and  the  channel  bed.  Bank 
resistance  per  unit  length  of  the  ice  cover  is  expressed  as  'zji  +  p/’j,  where  is  the  contribution  due 
to  cohesion  and  pf ^  is  the  ice-over-ice  friction  term,  with  being  the  stream  wise  force  per  unit  width 
of  the  cover.  Equation  52  was  derived  by  assuming  that  the  flow  conditions  when  the  ice  cover  reaches 
equilibrium  thickness  are  known.  Since  this  is  not  the  case,  a  solution  technique  taking  into 
consideration  the  interaction  between  the  flow  conditions  and  the  ice  cover  formation  was  developed 
(Lai  1988). 

Consider  the  case  of  an  equilibrium  ice  jam  condition  in  which  the  variables  do  not  change  along 
the  river.  Balancing  the  external  forces  acting  on  the  ice  cover  by  the  bank  resistance  yields 

2  (tcA  +  pi/)  Ax  =  (ti  +  Tg)  BAx  (53) 

where  /  =  longitudinal  stress  of  the  ice  cover 
Xj  =  bottom  friction  due  to  the  flow 

Xg  =  component  of  the  weight  in  the  direction  of  the  slope  of  the  cover 
p,  =  friction  coefficient  of  ice  to  the  banks 
Xj.  =  cohesive  strength  of  the  bank. 

Assuming  that  the  maximum  longitudinal  force  will  occur  at  a  passive  state,  it  can  be  expressed  as 
/■-Pil'-^)-^*!  ‘54) 

where  K2  =  ( 1  -  e)  tan^  ®  j 

<I>  =  internal  friction  angle  of  the  ice  accumulation 
e  =  porosity  of  the  ice  cover. 

The  friction  on  the  undersurface  of  the  ice  cover  can  be  expressed  as 

=  (55) 

and  the  weight  component  in  the  direction  of  the  slope  of  the  cover  can  be  expressed  as 

Xg  =  Pig/j5f.  (56) 

Using  eq  53-55  and  substituting  in  the  equilibrium  condition  of  eq  52,  we  obtain  the  following: 

F3(A,/i)=p^(i-^)/,2b2+2x^ 

^2  £i  =0  ^5.7^ 

^  7/3  [\  nj  p  A  \ 
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Start 


Figure  8.  Block  diagram  for  determining  the  initial  ice  cover  thickness. 


This  equation  has  to  be  solved  simultaneously  with  the  backwater  equation  (eq  51)  using  the 
Newton-Raphson  procedure  to  obtain  A  and  h.  The  starting  values  of  A  and  h  for  the  iterative  process 
can  be  obtained  either  from  the  previous  time  step  or  by  the  approximate  solution  obtained  from  eq 
52.  As  a  result  of  interaction,  progression  is  possible  over  a  wider  range  of  conditions  than  allowed 
by  eq  52.  A  simple  flow  chart  outlining  the  calculation  of  the  initial  ice  c  over  thickness  is  presented 
in  Figure  8. 

The  ice  cover  can  progress  as  a  result  of  the  obstruction  of  the  ice  Fow  due  to  natural  or  artificial 
causes.  At  an  existing  leading  edge  or  obstruction,  the  surface  fraction  of  the  incoming  ice  discharge 
will  collect  on  its  upstream  side  to  lengthen  the  cover.  The  computational  procedure  to  determine  the 
increase  in  ice  cover  length  can  be  illustrated  by  Figure  9.  Figure  9a  describes  the  condition  at  time 

r”.  The  concentration  profile  C{x,t'')  is  represented  by  A  ^A.  Figure  9b  describes  the  ice  condition  at 
During  the  time  interval  At,  the  water  parcel  A,  originally  located  at  the  leading  edge  X  l,  will 
move  to  A'.  The  curve  A '|A  represents  the  concentration  profile  at  r"^',  ignoring  production  and 
surface  accumulation.  This  profile  is  a  result  of  advection  of  the  profile  A ,  A.  Due  to  production  in  the 
open  water  area,  the  concentration  profile  at ignoring  surface  accumulation  only,  is 

represented  hyA'B  'C  ^2-  The  dotted  area  represents  the  ice  produced  during  At.  Since  an  fraction 
of  the  ice  that  is  passing  through  the  leading  edge  will  contribute  to  the  ice  cover  progression,  a  shaded 
&K&  A'C'D'E'  '\%  used  to  represent  the  ice  consumption  for  cover  progression.  In  the  Lagrangian 
scheme  the  profile  C  Y-V-t”'*'’ )  is  first  obtained  by  ignoring  the  consumption  of  surface  ice.  The  ice  in 
the  shaded  area  A'C'D'E'is  then  retrieved  to  form  the  ice  cover  between  X  i”  and  X  .  The  final 
concentration  profile  C(x,t''*^)  at  is  represented  by  A-^T)'E'A'. 
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a.  Condition  at  time  t". 


b.  Condition  at  time  l"*'. 

Figure  9.  Definition  sketch  for  ice  cover  progression. 


Letting  represent  the  surface  fraction  of  the  ice  passing  through  X  l  ,  the  length  of  the  progression 

can  be  obtained  by  equating  the  volume  of  the  ice  supply  for  the  cover  progression  +  a^^pCpAp 
and  the  volume  of  the  ice  in  the  cover  B/Up  (l-e^): 


X 


p 


_ Vc 

B/i(l  -Cc)  -  acCp/4p 


(58) 


where  h  =  thickness  of  the  newly  formed  ice  cover 
e^  =  porosity  of  the  ice  cover 
Cp  =  ice  concentration 
Ap  =  average  flow  area 

=  fraction  of  ice  discharge  going  into  ice  cover  formation. 


Ice  particles  that  remain  in  suspension  will  travel  under  the  ice  cover  after  they  pass  the  leading 
edge.  Because  of  the  insulation  of  the  ice  cover,  further  frazil  production  stops.  The  remaining  ice 
continues  to  rise  towards  the  underside  of  the  cover  under  the  influence  of  buoyancy  and  turbulent 
diffusion.  As  a  result,  ice  particles  will  be  deposited  along  the  underside  of  the  ice  cover  when  possible. 
Frazil  deposition  changes  the  ice  cover  thickness  and  influences  hydraulics.  During  the  deposition 
process,  the  flow  velocity  is  increasing  due  to  the  reduction  in  flow  area.  The  deposition  will  cease 
when  the  velocity  reaches  a  critical  velocity  of  deposition.  Deposited  frazil  ice  erodes  when  the  local 
flow  velocity  increases  beyond  a  critical  velocity  of  erosion.  Erosion  and  deposition  can  take  place  in 
different  parts  of  a  river  at  the  same  lime. 
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The  rate  at  which  the  frazil  ice  is  supplied  to  the  underside  accumulation  of  the  ice  cover,  due  to 
the  upward  movement  of  ice  in  the  suspension,  can  be  formulated  by  considering  the  mass  balance  of 
suspended  frazil  ice  mass: 

!(?£). -c.-.Sc‘'8  (59) 

dt  dx  ^ 

where  a‘^  is  the  probability  of  deposition  of  an  ice  particle  that  reaches  the  ice/water  interface  and 
is  the  buoyant  velocity  and  concentration  at  the  ice/water  interface.  Assuming  a (,C  =  and 
applying  the  continuity  equation,  eq  59  can  be  written  as 

=-avbCB.  (60) 

dt  dx 

In  Lagrangian  form  this  equation  becomes 

QC_  =_a  .  (61) 

Dt  A 

Hence,  the  decrease  in  ice  concentration  AC  in  the  suspension  over  a  travel  distance  Ax  is 


AC  =  Cl  1  -  exp  - 


av'bfi  A  x 


Equation  62  calculates  the  supply  of  frazil 
ice  to  the  underside  of  the  cover.  Deposition 
under  any  section  of  the  ice  cover  is  limited 
by  a  critical  flow  velocity  v  When  the 
thickness  of  the  frazil  accumulation  changes, 
the  hydraulic  condition  will  change  corre¬ 
spondingly  .  The  current  version  of  the  model 
assumes  that  the  water  level  of  the  river  does 
not  change  due  to  deposition  and  erosion,  but 
the  ice  cover  will  float  freely. 

Based  on  this  assumption  the  limiting 
condition  for  deposition  thickness  becomes 

hf<  (63) 

Pw\B  BVdepI 

where  A  =  flow  area 

hf  =  thickness  of  frazil  accumulation 
Pj^  =  density  of  frazil  accumulation 
o  =  subscript  representing  the  con¬ 
dition  at  the  beginning  of  the 
time  step. 


T  =  0‘’C 

(-)  I  (+) 


Erosion  of  frazil  ice  takes  place  when  the 
local  flow  velocity  over  a  frazil  deposition 
increases  beyond  a  certain  value  The 
limiting  thickness  of  frazil  ice  after  erosion 
can  be  determined  by  using  an  expression 
similar  toeq63: 


Figure  10.  Definition  sketch  for  ice  cover  growth. 
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(64) 


/if<  h(o  +  p-\^  --Q 
Pw  \  B  BVf 


Values  of  and  have  to  be  given  as  an  input  condition.  When  the  water  temperature  is  above 
freezing,  the  eroded  frazil  ice  will  first  be  melted  by  the  heat  of  the  water  before  bringing  down  the 
water  temperature. 

As  a  result  of  heat  exchanges  at  the  top  and  bottom  surfaces,  the  ice  cover  will  grow  or  decay  during 
the  winter.  A  finite-difference  model  developed  by  Shen  and  Lai  ( 1 986)  is  used  in  this  study.  In  this 
model  the  following  equation  is  used  to  calculate  the  rate  of  black  ice  growth  at  the  bottom  surface: 


PiL 


dhj 

dt 


Tm) 


where  and  h- 
and 

L 


=  thicknesses  of  snow  and  solid  ice 
=  thermal  conductivities  of  snow  and  ice 
=  latent  heat  of  fusion  for  ice  (Fig.  10). 


(65) 


If  frazil  ice  is  present  underneath  the  solid  ice  cover,  solid  ice  growth  will  occur  into  the  frazil  ice  layer, 
unaffected  by  turbulent  heat  exchange  from  the  river  flow.  The  rate  of  growth  of  solid  ice,  which  is 
higher  than  without  frazil  ice  accumulation,  can  be  expressed  as 

^  ,  dh.  T m  -  T, 

efPiL  — = - S! - 3- —  (66) 

dt  J_  +  ^  ^ 

^  wa  ^  i  ^  s 

where  is  the  porosity  of  the  frazil  ice.  Similarly  the  rate  of  change  of  the  frazil  ice  thickness  can  be 
expressed  as 


Pi£  ^  =  -  (7w  -  Tm)  -  etP.L  ^  .  (67) 

dt  ^  dt 

The  ice  growth  rate  in  the  absence  of  the  snow  layer  can  be  determined  by  letting  equal  0. 

Mechanical  failure  of  an  ice  cover  can  take  place  at  any  time  after  the  formation  when  the  internal 
strength  is  not  capable  of  withstanding  the  external  forces.  The  condition  for  ice  cover  failure  is 

2  (Xc/ +  Pi/)  <  (Tj  +  Xg -I- ta)B  (68) 

where/=  oi  +  o*  +  o-x,  and  oi ,  and  are  the  longitudinal  forces  per  unit  width  due  to  the  solid 
ice,  fragmented  ice  and  frazil  ice  layers,  respectively.  When  a  section  of  ice  cover  fails,  the  fragmented 
ice  masses  will  be  transponed  downstream,  where  they  will  accumulate  into  a  new  ice  cover  when 
conditions  permit. 


SIMULATION  OF  OHIO  RIVER  ICE  CONDITIONS 
Study  reach 

The  section  of  Ohio  River  simulated  in  this  study  extends  from  Pittsburgh,  Pennsylvania  (RM 
98 1.80),  to  Meldahl  Lock  and  D.'m,Ohio(RM  545.40),  with  a  total  distance  of  436.40  miles,  as  shown 
in  Figure  1 1 .  In  the  numerical  rr.odel  the  river  system  is  schematized  into  segments  connected  at  1 86 
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42' 


86'’  84‘>  82"  80'’  78” 


Figure  1 1 .  Ohio  River  between  Pittsburgh  and  Meldahl  L&D. 
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Belleville  (70) 
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Figure  12.  Schcmatization  of  the  upper  Ohio  River  system. 
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Table  1.  Schematizatiun  of  the  Ohio  River  system. 


Branch 

no. 

Upstream 

nodal 

point 

Downstream 

nodal 

point 

Description 

1 

1 

63 

Ohio  R.  from  Pittsburgh  to  Muskingum  R. 

2 

64 

85 

Ohio  R  from  Muskingum  R,  to  Kanawha  R. 

3 

86 

100 

Ohio  R.  from  Kanawha  R.  to  Big  Sandy  R. 

4 

tot 

111 

Ohio  R.  from  Big  Sandy  R.  to  Scioto  R. 

5 

1 12 

134 

Ohio  R.  from  Scioto  R.  to  Meldahl  L&D 

6 

l.3,S 

148 

Muskingum  R. 

7 

149 

162 

Kanawha  R. 

8 

163 

171 

Big  Sandy  R. 

9 

172 

186 

Scioto  R. 

Table 3.  Locks  and  damson  the  Ohio  River 
system. 


Pool  elevation 

Upstream 

Lock  and  dam 

(ft) 

nodal  point 

Table  2.  Gauged  lateral  inflow 

Emsworth 

Oashields 

710,0 

692.0 

4 

9 

distribution. 

Montgomery 

682.0 

17 

New  Cumberland 

664.5 

25 

Nodal  point 

Description  of  inflow 

Pike  Island 

644.0 

33 

Hannibal 

623.0 

46 

13 

Beaver  River 

Willow  Island 

602.0 

58 

66 

Little  Kanawha  River 

Belleville 

.582.0 

70 

69 

Hocking  River 

Racine 

560.0 

78 

96 

Guyandol  River 

Gallipolis 

538.0 

89 

99 

Twelve  Pole  Creek 

Greenup 

515.0 

108 

106 

Little  Sandy  River 

Meldahl 

485.0 

133 

1 1 1 

Little  Scioto  River 

Winfield 

566.0 

154 

nodal  points,  as  shown  in  Figure  12,  Some  of  the  tributaries  are  treated  as  separate  branches  (Table 
1 ),  and  others  are  treated  as  lateral  inflows  (Table  2).  The  entire  system,  which  covers  more  than  600 
miles,  consists  of  9  branches,  4  junctions  and  1 3  navigation  locks  and  dams  (Table  3 ).  The  nodal  points 
are  irregularly  spaced,  with  distances  ranging  from  1  to  5  miles,  with  closer  spaces  around  locks  and 
dams  or  junctions.  Figure  1 3  shows  the  discharges  of  the  Ohio  River  and  its  major  tributaries  for  the 
1985-86  winter,  provided  by  the  Ohio  River  Division  of  the  U.S.  Army  Corps  of  Engineers.  The 
geometric  data  tables  required  by  the  computer  model  RICEOH  were  constructed  using  data  furnished 
by  the  Ohio  River  Division  of  the  U.S.  Army  Corps  of  Engineers,  The  procedure  for  obtaining  the 
geometrical  data  was  fully  explained  by  Johnson  (1982). 

Lateral  inflows  were  furnished  by  the  Ohio  River  Division  as  either  gauged  or  ungauged  data.  The 
gauged  flows  represent  flow  from  small  rivers  that  are  not  treated  as  routed  branches.  The  ungauged 
data  are  normally  given  as  inflow  to  be  distributed  along  relatively  large  reaches  of  the  Ohio  River. 
Based  upon  navigation  maps  showing  natural  drainage  lines  into  the  Ohio  River,  the  distribution  of 
ungauged  lateral  inflows  presented  in  Table  4  has  been  prescribed  by  Johnson  ( 1982). 

Water  and  air  temperature  data  are  available  at  five  locks  and  dams:  Emsworth,  Montgomery, 
Hannibal,  Racine  and  Meldahl  (Fig.  14).  The  air  temperature  data  are  linearly  interpolated  to  each  node 
in  the  system.  The  water  temperature  at  the  upstream  boundary  on  the  Ohio  River  is  set  to  be  ()°C  when 
the  air  temperature  is  below  freezing  and  water  temperature  measurements  at  Emsworth  1.& Dare  close 
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a.  Discharges  of  the  Ohio  River  at  Node  No.  63  and  the  b.  Discharges  of  the  Ohio  River  at  Node  No.  88  and  the 
Muskingham  River  at  Node  No.  148.  Kanawha  River  at  Node  No.  162. 


c.  Discharges  of  the  Ohio  River  at  Node  No.  100  and  the  Big 
Sandy  Ri  ver  at  Node  No.  171 . 


d.  Discharges  of  the  Ohio  River  at  Node  No.  Ill  and  the 
Scioto  River  at  Node  No.  186. 


Figure  13.  Discharges  of  the  Ohio  River  and  its  major  tributaries  for  the  1985-86  winter. 


to  freezing;  otherwise  it  was  assumed  to  be  the  same  as  the  water  temperature  at  Ems  worth  L&D.  Since 
no  data  were  available  for  tributaries,  the  air  temperature  along  the  tributary  is  assumed  to  be  the  same 
as  the  air  temperature  at  the  junction  of  the  tributary  and  the  main  stream.  The  water  temperature  at 
the  upstream  end  of  each  tributa.y  was  set  to  be  the  same  as  the  water  temperature  immediately 
upstream  of  the  junction  of  the  tributary  and  the  main  stream. 


Numerical  simulation  and  results 

The  computer  model  was  appl  ied  to  the  upper  Ohio  River  system  for  the  1 985-86  winter.  The  first 
step  in  the  simulation  was  to  calibrate  the  heat  exchange  coefficient  h^^.  The  calibration  was  made 
based  on  the  data  at  Montgomery  and  Hannibal  using  a  least-squares  technique  (Lai  1988).  These  two 
stations  were  selected  since  there  is  no  tributary  inflow  between  them.  This  calibration  gives  h^^=23.9 
Wm~^°C“’ .  This  value  is  applied  to  the  entire  study  reach.  The  field  data  also  indicate  that  the  water 
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Temperature  ('C) 


Table  4.  Ungauged  lateral  inflow  distribution. 

Nodal  Percent  Nodal  Percent 


point 

Unhanged  reach 

of flow 

point 

Ungaiiged  reach 

of flow 

3 

Pitlsburgh-Wheeling 

5 

81 

Pomeroy-Huntington 

20 

8 

Pillsburgh-Whceling 

5 

83 

Pomeroy-Huntington 

10 

10 

Pitlsburgh-Wheeling 

5 

84 

Pomeroy-Huntington 

10 

13 

Pitlsburgh-Wheeling 

5 

87 

Pomeroy-Huntington 

10 

16 

Pitlsburgh-Wheeling 

7 

91 

Pomeroy-Huntington 

20 

20 

Pitlsburgh-Wheeling 

13 

97 

Pomeroy-Huntington 

30 

24 

Pitlsburgh-Wheeling 

13 

100 

28 

Pitlsburgh-Wheeling 

17 

105 

Huntington-Maysville 

30 

32 

Pitlsburgh-Wheeling 

30 

114 

Huntington-Maysville 

10 

100 

120 

Huntington-Maysville 

30 

36 

Wheeling-St.  Mary 

18 

124 

Huntington-Maysville 

30 

40 

Wheeling-St.  Maty 

13 

100 

45 

Wheeling-St.  Mary 

16 

129 

Maysville-Cincinnati 

30 

51 

Wheeling-St.  Mary 

30 

132 

Maysville-Cincinnati 

40 

55 

Wheeling-St.  Mary 

23 

N/A 

Maysville-Cincinnati 

30 

100 

100 

67 

St.  Mary-Pomeroy 

10 

68 

St.  Mary-Pomeroy 

10 

73 

St.  Mary-Pomeroy 

30 

76 

St.  Mary-Pomeroy 

20 

77 

St.  Mary-Pomeroy 

10 

77 

St.  Mary-Pomeroy 

20 

1  "'O 

a.  Emsworth  L&D.  h.  Montgomery  L&D. 

Figure  14.  Air  and  water  temperatures  at  five  locks  and  dams. 
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25 


c.  Hannibal  L&D. 


d.  Racine  L&D. 


e.  Meldahl  L&D. 

Figure  14  (coni' d).  Air  and  water  temperatures  at  five  locks  and  dams. 
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Figure  15.  Comparison  of  the  water  temperatures  at  Emsworth  L&D, 
Montgomery  L&D  and  the  Beaver  River. 


temperature  at  Montgomery  is  approximately  2°C  higherthan  that  of  Emsworth,  although  Montgom¬ 
ery  is  only  25  miles  downstream  of  Emsworth  and  the  air  temperatures  at  these  two  stations  are  about 
the  same.  This  phenomenon  can  only  be  explained  by  considering  thermal  effluents  discharged  into 
the  river  by  local  industries.  Since  there  is  no  record  of  the  magnitude  and  distribution  of  thermal 
effluents,  it  is  assumed  that  the  effect  of  thermal  effluents  can  be  modeled  by  considering  thermal 
effluent  discharges  into  the  Ohio  River  from  the  Beaver  River.  The  validity  of  this  assumption  is 
substantiated  by  the  existence  of  a  large  open-water  area  in  the  vicinity  of  the  Beaver  River  and  Ohio 
River  confluence  during  the  midwinter  ice-covered  period  (Gatto  et  al.  1987).  A  calibration  was  then 
carried  out  to  determine  the  water  temperature  of  the  Beaver  River  inflow.  Figure  15  shows  the 
calculated  Beaver  River  water  temperature  and  the  measured  and  calibrated  water  temperatures  at 
Montgomery  and  Emsworth.  The  temperature  of  the  Beaver  River  inflow  is  approximately  4°C  higher 
than  the  Ohio  River  temperature  during  the  calibration  period.  This 
difference  is  expected  to  decrease  during  the  winter.  It  is  recom¬ 
mended  that  additional  water  temperature  measurements  be  obtained 
at  the  downstream  end  of  the  Beaver  River  to  be  able  to  accurately 
account  for  the  effect  of  thermal  effluents  in  this  reach.  Figure  16 
shows  the  comparison  between  the  observed  and  simulated  water 
temperatures  at  Hannibal,  along  with  the  observed  water  temperature 
at  Emsworth.  Similar  comparisons  for  water  temperatures  at  Racine 
and  Meldahl  L&D  are  shown  in  Figure  17. 

Tosimulate  ice  cover  conditions,  several  additional  parameters  are 
needed.  These  include  the  ice  floe  thickness,  the  underside  roughness 
of  the  ice  cover,  the  critical  Froude  numbers,  the  buoyant  velocity  and 
the  fraction  of  the  surface  ice  discharge  that  will  be  used  to  form  the 
ice  cover.  The  values  of  these  parameters  used  in  the  simulation  are 
summarized  in  Table  5.  The  upstream  boundary  discharges  of  the 
Ohio  River  and  its  tributaries  during  the  simulation  period  are  shown 


Table  5.  Parameter  values 
for  sample  simulation. 


Parameter 

Value 

'‘wa 

1 

23.9  W  m--  °C-' 

0.08  m 

", 

0.015  or  0.020 

0.06 

F,  , 

0.09 

r,max 

^'dep 

0.6  m/s 

'cro 

0.7  m/s 

“''b 

0.001  m/s 

“c 

0.85 

0.2 

‘V 

0.6 

u 

0.28 

0.98  kPa 
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November  1985 


Figure  16.  Comparison  of  the  water  temperatures  at  EmsworthL&D  and 
Hannibal  L&D. 


November  1985 

Figure  17.  Comparison  of  the  water  temperatures  at  Emsworth  L&D, 

Racine  L&D  and  Meldahl  L&D. 

in  Figure  1 8.  The  rating  curve  used  at  the  downstream  boundary  (Node  No.  1 34)  was  constructed  from 
the  data  supplied  by  the  Ohio  River  Division  (Fig.  19).  The  underside  roughness  of  an  ice  cover  can 
affect  the  energy  slope  and  therefore  the  water  levels.  Since  no  observed  water  level  data  are  available 
for  the  Ohio  River,  no  calibration  was  made  and  the  roughness  coefficient  is  assumed  to  be  equal 
to  0.015  or  0.02. 

The  ice  cover  on  the  Ohio  River  is  formed  only  by  juxtaposition,  and  no  narrow  or  wide  jams  were 
observed  either  in  the  field  or  in  the  simulated  result.s.  The  thickness  of  the  initial  ice  cover  is  the  main 
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Days  from  Dec  21 .  1985 


Figure  18.  Discharges  during  the  simulation  period. 


Discharge  (ll^/s) 

Figure  19.  Rating  curve  at  the  downstream  boundary. 

parameter  governing  the  rate  of  ice  cover  progression.  In  general,  the  ice  cover  progresses  upstream 
very  fast  immediately  after  its  initiation.  During  warm  periods  the  cover  can  melt  away  completely 
without  re-forming. 

Two  sets  of  data  were  available  for  comparison  between  the  simulated  and  observed  ice  cover 
lengths.  The  ice  atlas  for  the  Ohio  River  1985-86  .sea.son  (Gatto  et  al.  1987)  has  a  few  days  of 
observations  during  the  winter  when  the  river  was  covered  with  ice.  The  ice  atlas  was  prepared  from 
videotapie  records.  In  addition,  navigation  charts  based  on  manual  estimates  of  the  lengths  of  ice  covers 
from  lock  and  dams  were  available.  The.se  data  are  less  reliable  and  do  not  compare  well  with  the  ice 
atlas.  The  navigation  chart  data  becomes  inaccurate  when  the  length  of  the  ice  cover  was  greater  than 
the  visible  distance  from  the  dam  and  smaller  than  the  length  of  the  pool.  However,  these  charts  provide 
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a.  Ice  cover  roughness  =  0.015  and  one  thermal  lateral  inflow. 
Figure  21 .  Observed  and  simulated  ice  covers. 


good  information  on  the  first  appearance  of  the  ice  cover  in  each  pool,  and  this  information  compares 
well  with  the  simulation,  as  shown  in  Table  6. 

Preliminary  calibration  runs  show  that  there  is  a  significant  heat  inflow  between  Montgomery  and 
New  Cumberland  L&D.  This  influx  of  heat  affects  the  length  of  the  ice  cover  in  the  New  Cumberland 
pool.  In  the  simulation  this  heat  influx  is  accounted  for  by  increasing  the  water  temperature  of  lateral 
inflows,  as  shown  in  Figure  20.  Figure  21  shows  the  comparison  of  simulated  and  observed  lengths 
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Emswotih 

Dashields 

Montgomery 

New  Cumberland 
Pike  Island 

Hannibal 
Willow  Island 

Belleville 

Racine 


Galllpolis 


Greenup 


Meldahl 


b.  Ice  cover  roughness  =  0.015  and  two  thermal  lateral  inflows. 


Figure  21  (cont'd). 
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Table  6.  First  appearance  of  the  ice  cover. 


of  ice  covers  in  pools  along  the  Ohio  River.  In  Figure  21a 
a  water  temperature  of  4°C  was  assigned  to  one  of  the 
lateral  inflows  between  Montgomery  and  New  Cumberland 
L&D.  In  Figures  21b  and  21c  the  4°C  water  temperature 
was  assigned  to  both  lateral  inflows.  Figures  21b  and  21c 
show  that  the  effect  of  is  relatively  insignificant. 


SUMMARY  AND  CONCLUSION 


Lock  and  dam 

Simulated 

Ohsened 

Emswonh 

nrit 

\iai 

Dashields 

12/26 

12/27 

Monlgomer)' 

12/27 

12/27 

New  Cumberland 

12/27 

12/26 

Pike  Island 

12/27 

12/26 

Racine 

12/28 

12/27 

The  purpose  of  this  study  was  to  develop  a  model  for  simulating  river  ice  and  flow  conditions.  The 
model  ORICE  was  developed  and  tested  on  the  Ohio  River  system  between  Pittsburgh  and  Meldahl. 
Together  with  good  weather  forecasts  and  proper  geometric  data  required  by  program,  this  model  can 
be  applied  to  any  river  system  of  dendritic  configuration  to  forecast  river  ice  conditions. 

Field  data  on  the  Ohio  River  system  are  marginally  adequate.  Comparisons  between  the  lengths  of 
simulated  and  obser\'ed  ice  covers  are  given  in  Figure  21.  We  can  conclude  that  the  simulation 
compares  very  well  with  the  ice  atlas. 

The  information  on  the  ice  cover  conditions  from  navigation  charts  and  simulations  compare  very 
well  in  pools  of  the  lower  dams  and  those  of  Emsworth  and  Dashields.  Slight  discrepancies  can  be 
noticed  in  the  Montgomery  and  New  Cumberland  pools.  Significant  differences  between  the  ice  atlas 
and  the  navigation  charts  are  also  found  in  these  two  reaches.  This  poses  a  reasonable  doubt  as  to  the 
accuracy  of  the  information  from  the  navigation  charts.  The  accuracy  of  the  simulation  can  be 
improved  if  more  data  on  thermal  discharges  can  be  obtained.  We  believe  that  the  thermal  discharges 
that  were  assumed  in  this  study  are  not  the  only  ones  on  the  Ohio  River.  Also,  correct  water 
temperatures  at  the  upstream  boundaries  of  tributaries,  together  with  air  temperature  measurements, 
would  provide  information  for  more  accurate  simulation. 

To  calibrate  the  ice  cover  roughness  accurately,  it  would  be  necessary  to  collect  water  stage  data 
and  compare  these  with  the  simulation.  The  present  study,  however,  shows  that  ice  cover  progression 
is  not  very  sensitive  to  variations  in  ice  cover  roughness  /tj. 
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APPENDIX  A:  COMPUTER  PROGRAM  AND  USER’S  MANUAL 


The  model  RICEOH  consists  of  two  major  pans;  the  unsteady  flow  computation  and  the  water 
temperature  and  ice  cover  condition  simulations.  The  unsteady  flow  computation  is  based  on  the  model 
set  up  and  tested  by  Johnson  ( 1 982).  Minor  changes  are  made  to  account  for  the  ice  cover  effect.  The 
water  temperature  and  ice  simulation  model  was  based  on  the  river  ice  model  developed  by  Lai  ( 1 988). 
The  model  of  Lai  was  originally  developed  for  a  single-channel  river.  Improvements  are  made  in  the 
present  model  to  accommodate  a  dendritic  system.  In  addition  RICEOH  can  treat  lateral  thermal 
inflows.  The  program  was  written  in  Fortran??.  Different  ice  processes  are  modeled  separately  in  the 
form  of  independent  modules  that  are  convenient  for  future  modifications.  The  flow  chart,  which 
outlines  the  computational  steps,  is  given  in  Figure  Al. 
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Figure  Al.  Flow  chart  of  the  computational  steps  in  RICEOH. 
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Subroutine  descriptions 

The  program  is  composed  of  a  main  program  that  controls  the  flow  of  computations  and  30 
subroutines  that  perform  various  functions.  A  brief  description  of  the  role  of  each  subroutine  is 
presented  below.  All  of  the  subroutines  from  the  original  unsteady  flow  program  (Johnson  1982) 
remain  essentially  the  same. 

LOCKDA  sets  coefficients  in  the  equations  applied  to  a  computational  reach  containing  a 
navigation  lock  and  dam  so  that  and  =  /{r),  where  /is  a  known  function  of  time. 

LINEAR  linearly  interpolates  values  of  the  flow  area,  top  width  and  Manning’s  roughness 
coefficient  for  a  particular  water  surface  elevation  from  geometric  tables. 

FLOOD  determines  the  surface  area  of  the  flood  plain  used  in  computing  the  exchanges  between 
the  main  channel  and  the  flood  plain  (^|)|. 

INITIA  initializes  the  hydraulic  computations.  Various  minor  computations  are  performed,  e.g. 
the  flood  plain  area,  which  is  input  as  a  cross-sectional  area  and  then  changed  to  surface  area,  and  the 
spatial  computational  steps  are  computed. 

JOINTF  joins  the  results  from  the  forward  sweep  on  the  main  river  and  a  tributary  with  the 
confluence  equations  so  that  the  forward  sweep  on  the  main  river  can  continue. 

COEFFI  computes  the  coefficients  in  the  system  of  linear  algebraic  equations  (eq  24  and  25), 
which  depend  on  information  known  from  the  previous  time  step.  These  coefficients  are  computed  in 
this  subroutine. 

FORWAR  computes  the  coefficients  in  the  forward  sweep  of  the  double-sweep  algorithm. 

BACWAR  computes  the  unknowns  in  the  backward  sweep  of  the  double-sweep  algorithm. 

NEWFLO  is  called  at  the  end  of  each  time  step  to  update  flow  conditions  for  initiating 
computations  in  the  next  time  step. 

DAMRC  computes  pool  elevations  upstream  of  the  lock  and  dam  from  a  rating  curve  at  the 
structure. 

DOWNCO  computes  difference  coefficients.  The  rating  curve  at  the  downstream  end  of  the  main 
river  is  given  in  the  form  of  a  number  of  linear  segments.  Each  linear  segment  is  defined  by  specifying 
the  discharge  and  corresponding  elevation  at  the  end  of  the  segment.  DOWNCO  uses  this  information 
to  compute  the  difference  coefficients. 

BRYC  AL  inputs  time-dependent  boundaiy  data  and  returns  either  linearly  interpolated  data  or  the 
data  that  have  just  been  read. 

LATERA  performs  the  same  function  as  BRYCAL  except  that  lateral  inflows  instead  of  open¬ 
boundary  data  are  returned. 

OPEN  defines  all  input  and  output  files. 

REINT  reads  all  information  needed  for  thermal  computations. 

LAL  links  the  main  program  and  thermal  computations. 

PREPAR  creates  a  separate  subsystem  for  each  branch  of  the  river  system.  All  computations  will 
be  performed  on  this  subsystem. 

POVRAT  returns  results  of  thermal  computations  from  the  subsystem  to  the  main  system. 

ICE  performs  all  ice-related  computations  by  calling  the  subroutines  for  individual  ice  processes 
such  as  progression,  erosion  and  growth. 

AVERAG  computes  the  reach-average  values  of  geometrical  and  hydraulic  variables  from  the 
known  nodal  values.  TTiis  is  required  because  all  hydraulic  computations  are  performed  using 
variables  defined  at  nodal  points,  while  ice  routines  require  reach-averaged  variables. 

NITHI  computes  the  initial  thickness  of  the  ice  cover  for  a  given  reach  when  the  reach  number  is 
provided  as  a  dummy  argument.  This  subroutine  has  to  be  called  once  for  each  reach  that  has  a  leading 
edge  and  for  which  an  initial  ice  Ihickne.ss  has  to  be  computed.  The  initial  thickness  mainly  depends 
on  the  hydraulics  conditions  of  the  reach,  and  in  the  process  of  detennining  it.  juxtaposition,  narrow 


34 


and  wide  jam  modes  are  considered.  Ice  cover  progression  is  limited  by  the  maximum  critical  Froude 
number. 

LAGRAN  computes  the  distribution  of  water  temperature  or  frazil  ice  concentration  along  the 
river.  It  uses  a  one-dimensional  dispersion  equation  after  neglecting  the  longitudinal  dispersion  term. 
A  Lagrangian  method  is  used  to  solve  the  equation.  This  method  assumes  that  parcels  of  water  having 
known  temperatures  are  released  from  nodal  points  at  the  beginning  of  every  time  step.  Their  locations 
and  concentrations  are  computed  at  the  end  of  every  time  step.  Nodal  temperatures  are  determined 
from  these  values  using  linear  interpolation. 

OPENLA  is  called  for  each  open  water  area.  The  results  are  stored  separately  to  be  used  for  the 
following  reasons:  a)  When  the  subroutine  LAGRAN  is  called,  the  final  frazil  concentration 
distribution  is  the  sum  of  the  convected  ice  and  the  ice  generated  in  the  open-water  patches;  b)  When 
the  volume  of  ice  used  for  deposition  or  ice  cover  progression  at  a  given  reach  is  computed,  frazil  ice 
suspension,  which  has  already  completed  the  motion,  is  pulled  back  to  add  to  the  volume.  OPENLA 
helps  to  sort  out  the  amounts  of  ice  created  in  different  open-water  reaches  to  make  sure  that  any 
suspended  portion  of  ice  is  deposited  upstream  from  where  it  was  produced. 

PROGRE  is  called  at  every  leading  edge  of  the  ice  cover  once  during  every  time  step.  The  current 
location  of  the  leading  edge  has  to  be  provided  as  a  dummy  argument.  The  subroutine  computes  the 
volume  of  ice  that  contributes  to  cover  formation  by  taking  a  fraction  of  the  ice  discharge  that  has 
passed  the  leading  edge.  PROGRE  computes  all  ice  cover  progressions  that  can  take  place  at  leading 
edges  during  the  time  step.  Both  the  ice  thickness  and  the  horizontal  extent  of  progression  are 
computed.  Changes  in  the  local  hydraulics  are  computed  based  on  a  backwater  approximation. 

OPLA  uses  the  total  amount  of  ice  computed  in  LAGRAN  to  determine  the  amount  of  ice  that  can 
contribute  to  deposition.  It  separates  the  amounts  of  ice  produced  at  open- water  reaches  downstream 
of  the  current  reach  because  this  ice  cannot  be  deposited  in  any  upstream  reach.  This  separation  is 
necessary  because  deposition  or  progression  is  computed  as  a  reverse  Lagrangian  process  or  “pulling 
back”  process. 

TRAPEZ  computes  nodal  values  of  ice  concentration.  During  the  computation  of  the  ice  cover 
progression  and  deposition,  a  portion  of  the  moving  ice  is  removed.  Since  concentrations  are  assigned 
to  the  nodes,  the  amount  of  ice  removed  from  a  reach  takes  a  shape  of  a  trapezoid  in  the  longitudinal 
concentration  profile.  After  removal  of  these  trapezoidal  shapes,  the  remaining  concentrations  have 
to  be  converted  to  nodal  values  before  beginning  of  the  next  time  step.  Subroutine  TRAPEZ  computes 
these  nodal  values  after  conserving  the  mass  and  convection.  This  subroutine  is  called  by  PROGRE 
and  ERODEP  subroutines. 

ERODE?  computes  the  amount  of  erosion  and  deposition  that  takes  place  in  a  given  reach.  It  is 
called  once  during  every  time  step  for  every  ice-covered  reach  in  the  river.  The  occurrence  of  erosion 
and  deposition  in  a  river  reach  is  decided  depending  on  the  velocity  of  flow.  In  the  computation  of 
deposition,  the  total  ice  quantity  that  can  be  deposited  is  computed  first.  After  depositing  up  to  a  critical 
limit,  the  excess  ice  is  put  back  in  the  stream  as  concentration.  Erosion  is  computed  similarly,  using 
a  critical  velocity  criterion.  Since  the  rate  of  erosion  is  not  known,  the  eroded  ice  is  distributed 
uniformly  over  the  travel  distance  of  a  particle.  If  the  water  is  warm,  part  or  all  of  the  ice  cover  can  be 
melted. 

SHOVE  checks  the  stability  of  the  ice  cover  against  the  net  streamwise  gravity  and  friction  forces. 
The  strength  of  the  solid  ice  crust  is  considered  in  the  computations.  Parts  of  the  initial  ice  cover  and 
frazil  deposited  are  considered  as  granular  materials  under  passive  stress  conditions.  If  external  forces 
are  greater  than  the  strength,  the  ice  cover  is  assumed  to  fail  in  the  current  reach,  as  well  as  the  upstream 
reaches  downstream  of  the  next  ice  control  structure  or  ice  bridge. 

PUTBAC  is  called  only  if  the  ice  cover  in  a  given  reach  is  known  to  fail  as  indicated  by  SHOVE. 
It  first  tries  to  re-form  an  initial  ice  cover  using  the  fragments  of  broken  cover.  The  formation 
conditions  and  the  thickness  of  this  initial  cover  are  decided  by  the  local  hydraulic  conditions.  After 
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shoving,  a  new  thicker  ice  cover  can  form,  or  the  broken  pieces  can  go  under  the  cover  and  add  to  the 
ice  discharge.  This  ice  discharge  travels  downstream  and  can  form  a  cover,  be  deposited  under  an 
existing  cover  or  flush  away  through  the  downstream  boundary.  Subroutine  PUTBAC  is  called  once 
during  every  time  step  for  every  reach  that  is  subjected  to  shoving. 

GROWTH  is  capable  of  simulating  the  thermal  growth  and  decay  of  an  ice  cover  in  a  given  river 
reach.  It  has  to  be  called  once  during  every  time  step  for  every  river  reach  that  is  partially  or  fully 
covered  with  ice.  The  model  uses  a  steady-state  finite-difference  method  to  compute  incremental 
growth  or  decay  during  every  time  step.  Heat  exchanges  at  the  top  and  bottom  surfaces  are  computed 
using  linearized  models.  The  model  can  simulate  both  black  and  white  ice  growth  while  having  a  layer 
of  snow  on  top  and  frazil  ice  at  the  bottom. 

Input  data 

Input  data  for  the  model  can  be  divided  into  two  major  groups,  i.e.  input  data  for  hydraulic 
computations  and  input  data  for  thermal-ice  computations.  Input  data  for  the  hydraulics  part  consist 
of  four  data  files.  FILE4.DAT  contains  general  information  about  the  system,  such  as  the  number  of 
nodal  points,  the  number  and  description  of  branches  and  joints  and  description  of  dams  and  lateral 
inflows.  FILE2.DAT  contains  prescribed  initial  stages  and  discharges  for  every  point  of  the  system 
and  discharges  for  dams  that  had  a  rating  curve  as  input  in  nLE4.DAT.  FILE5.DAT  contains  all 
geometrical  information,  such  as  cross-sectional  area,  top  width.  Manning’s  roughness  coefficient  and 
flood  plain  area,  given  as  a  function  of  stage  for  every  point  of  the  system.  FILE  1  .DAT  contains  daily 
values  of  lateral  inflows  for  every  lateral  inflow  point,  boundary  discharges  for  every  branch  in  the 
system,  and  elevations  for  all  dams  in  the  system  that  do  not  have  a  rating  curve  as  input.  These  files 
are  the  same  as  files  used  in  the  unsteady  flow  model  for  the  Ohio  River  (Johnson  1982)  except  that 
the  system  immediately  downstream  of  Meldahi  was  not  included. 

Input  data  for  thermal-ice  computations  are  separated  into  two  data  files.  FILE6.DAT  contains 
general  parameters  used  in  thermal  computations,  including  the  density  of  water  and  ice  or  the  water- 
to-air  heat  exchange  coefficient,  followed  by  general  ice  parameters  for  every  point,  namely,  ice 
roughness,  particle  and  critical  velocities  and  ice  bridging  description.  The  boundary  temperature 
conditions  are  also  stored  in  this  file.  During  each  day  of  the  simulation,  air  temperatures  at  points 
where  this  information  is  available  are  prescribed.  This  includes  temperatures  at  five  lock  and  dams 
and  temperatures  at  the  first  upstream  and  last  downstream  points  of  each  branch.  Also,  boundary 
water  temperatures  for  every  branch  and  lateral  inflow  that  is  important  for  thermal  computations  are 
stored  in  this  file. 

FILE3.DAT  contains  the  initial  temperatures  for  both  air  and  water  at  every  point  of  the  system. 
Depending  on  whether  an  ice  cover  already  exists  in  the  system  or  not,  it  may  provide  descriptions  of 
present  ice  conditions  for  every  river  reach.  This  includes  the  length  of  the  ice  cover,  the  thickness  of 
the  solid  ice  cover,  the  frazil  deposition  thickness,  the  submerged  thickness  and  the  initial  ice  cover 
thickness.  Water  temperatures  in  this  file  can  be  negative.  However,  this  means  that  there  is  a  frazil 
ice  concentration  at  that  point,  and  the  water  temperature  is  actually  0°C. 

Detailed  explanations  of  every  input  statement  in  the  program  and  a  sample  of  input  data  files 
follow.  More  information  about  input  data  for  the  hydraulic  part  of  the  program,  as  well  as  plotting 
capabilities,  is  provided  by  Johnson  (1982). 

Input  statements 

Statement  I 

Variable:  TITLE(I) 

Format:  10A8 
File  read:  FILE4 
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Variable 


Value 


Meaning 


TITLE  Descriptive  title  of  the  run 


Statement  2 

Variables:  IGEOM,  ILUG,  IPUNCH,  IBACK 
Format:  415 
File  read:  nLE4 


Variable 

Value 

Meaning 

IGEOM 

0 

Geometry  tables  printed 

1 

Geometry  tables  not  printed 

ILUG 

90 

File  from  which  geometry  tables  are  read 

IPUNCH 

vtO 

Initial  conditions  not  printed 

IBACK 

Less  printed  output 

Statement  3 

Variables:  NC,  NBRS,  NJUNC,  NDAMS,  NXMAIN,  1ST  AGE,  IFPLN,  TOTALT,  TST 

Format:  715,  lOX,  2F10.0 

File  read:  nLE4 

Variable 

Value 

Meaning 

NC 

185 

Total  number  of  nodal  points  minus  one 

NBRS 

9 

Total  number  of  branches 

NJUNC 

4 

Total  number  of  junctions 

NDAM 

13 

Total  number  of  dams 

NXMAIN 

134 

Last  nodal  point  on  the  main  stream 

ISTAGE 

5 

Number  of  entries  in  channel  geometry  tables 

IFPLN 

3 

Number  of  entries  in  flood  plain  tables 

TOTALT 

2 

Number  of  days  of  computation 

TSTEP 

3600 

Time  step  in  seconds 

Statement  4 

Variables:  NSTAT,  IPRINT,  INTVG,  INTVD,  INTVP,  NOXS 

Format:  615 

File  read:  FILE4 

Variable 

Value 

Meaning 

NSTAT 

37 

Number  of  nodal  points  at  which  output  is  requested 

IPRINT 

5t0 

Detailed  output  provided 

INTVG 

24 

Major  print  interval 

INTVD 

24 

Small  print  interval  for  particular  days  (see  following  card) 

INTVP 

24 

Interval  for  placing  points  in  plot  file 

NOXS 

0 

Number  of  stations  at  which  plots  are  desired 

37 


Statement  5 

Variables:  STOP,  ISDDP,  ISMDP,  ISYDP,  ETDP,  lEDDP,  lEMDP,  lEYDP 
Format:  FIO.O,  315, 5X,  FIO.O,  315 
File  read:  FILE4 


Variable 

Value 

Meaning 

STOP 

Starting  time  on  24-hr  clock  for  small  print  interval 

ISDDP 

Starting  day  for  small  print  interval 

ISMDP 

Starting  month  for  small  print  interval 

ISYDP 

Starting  year  for  small  print  interval 

ETDP 

Ending  time  on  24-hr  clock  for  small  print  interval 

lEDDP 

Ending  day  for  small  print  interval 

lEMDP 

Ending  month  for  small  print  interval 

lEYDP 

Ending  year  for  small  print  interval 

If  INTVD  =  INTVG,  irrsert  blank  card  for  card  5. 


Statement  6 

Variables:  NPRINT(I), 

I=1,NSTAT 

Format:  1615 

File  read:  FILE4 

Variable  Value 

Meaning 

NPRINTd) 

Nodal  points  at  which  output  is  desired 

Statement  7 


Variables:  IPLT,  ISCL,  lYCN,  lOPl 

Format:  415 

File  read:  FILE4 

Variable  Value  Meaning 

IPLT 

0 

No  plots 

1 

Elevation  plots 

2 

Discharge  plots 

3 

Elevation  and  discharge  plots 

4 

Velocity  plots 

5 

Elevation  and  velocity  plots 

6 

Discharge  and  velocity  plots 

7 

Elevation,  discharge  and  velocity  plots 

ISCL 

1440 

Number  of  minutes/iiKh  on  jc-axis 

lYLN 

0 

Length  of  y  axis  in  inches  (0  defaults  to  8  inches) 

lOPl 

0 

Different  y  interval  on  each  plot 

1 

Same  y  interval  on  each  plot 
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Statement  8 

Variables;  STITLE(I),  NPLOTa) 
Format:  A32, 3X,  15 
File  read:  FILE4 

Variable  Value  Meaning 


STITLE(I)  Title  of  station 

NPLOT(I)  Nodal  point  number  of  /***  station 


Statement  9 

Variables:  ID,  IBRNCHa.l),  IBRNCH(1, 2),  IBRS(I),  IBC(I) 
Format:  515 
File  read:  FILE4 

Variable  Value  Meaning 


ID 

IBRNCH(I.l) 
IBRNCH(I.2) 
IBRS(I)  1 

0 
-1 

IBC(I)  -1 

0 
1 


Branch  number 
First  nodal  point  of  the  branch 
Last  nodal  point  of  the  branch 
Branch  has  upstream  outer  boundary 
Interior  branch 

Branch  has  downstream  outer  boundary 

Rating  curve  will  be  used  if  this  branch  has  downstream  outer  boundary 
This  is  an  interior  branch 

Elevations  will  be  input  if  this  branch  has  downstream  outer  boundary 


There  is  one  card  no.  9  for  each  branch.  All  mainstream  branches  should  be  numbered  before 
numbering  tributaries. 


Statement  10 

Variables:  ID,  UUNC(I,K),  K=l,3,  AL(J),  AL(J+1) 
Format:  415,  lOX,  2F10.0 
File  read;  FILE4 

Variable  Value  Meaning 


ID 

UUNC(I,1) 

UUNC(I,2) 

UUNC(I,3) 

AL(J) 


Junction  number 

Branch  number  of  the  main  river  upstream  of  the  Junction  7 
Tributary  branch  number  upstream  of  the  junction  I 
Branch  number  of  the  main  river  downstream  of  the  Junction  I 
Energy  correction  coefficients  corresponding  to  node  7 


There  is  one  card  no.  10  for  each  Junction. 
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Statement  11 

Variables;  TDAM(I),  HSETO(I),  NL(I),  NVARY(I) 
Format:  A8, 2X,  FIO.O,  315 
File  read:  FILE4 


Variable  Value 

Meaning 

TDAM(I) 

HSETO(I) 

NL(I) 

NVARY(I)  0 

1 

2 

Name  of  the  dam 

Pool  elevation  maintained  by  dam 

Nodal  point  immediately  upstream  of  the  dam 
Pool  elevation  equal  to  HSETO(I) 
Time-varying  elevations  of  pool  are  input 
Rating  curve  is  input  for  P**  dam 

Statement  12 

Variables:  lEDYHD 

Format:  15 

File  read:  FILE4 

Variable  Value 

Meaning 

lEDYHD  0 

Eddy  head  loss  coefficients  set  to  0 

If  lEDYHD  =  0,  skip  the  following  card. 

Statement  13 

Variables:  N,  CKE(N) 
Format:  15,  FI 0.2 

File  read;  nLE4 

Variable  Value 

Meaning 

N 

CKE(N) 

Number  of  nodal  point 

Eddy  head  loss  coefficient  for  nodal  point  N 

Card  no.  12  is  repeated  lEDYHD  times. 


Statement  14 

Variables:  NRCH,  (IRCH(I),  1=1,  NRCH) 
Format:  1615 
File  read;  FILE4 

Variable  Value  Meaning 


NRCH  38  Total  number  of  reaches  containing  lateral  inflows 

IRCH(I)  Numbers  of  the  reaches  containing  lateral  inflow 
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Statement  15 

Variables:  MNTH.  KDAY,  K't^AR,  TIME 
Format:  315,  5X,  FIO.O 
File  read:  FILE2 

Variable  Value  Meaning 


MNTH 

12 

Starting  month 

KDAY 

21 

Starting  day 

KYEAR 

86 

Starting  year 

TIME 

7 

Starting  time  on  a  24-hr  clock 

Statement  16 

Variables:  H(I,JSP),  1=1,  NX 
Format:  8F10.0 
File  read:  FILE2 

Variable  Value  Meaning 


H(IJSP) _ Initial  water  surface  elevation  in  feet  at  each  nodal  point 


Statement  17 

Variables:  Q(I,JSP),  1=1,  NX 
Format:  8F10.0 
File  read:  FILE2 

Variable _ Value  Meaning 


Q(I,JSP) 


Initial  discharge  in  cfs  at  each  nodal  point 


Statement  18 

Variables:  QL2L(K),  K=l,  NRCH 
Format:  8F10.0 
File  read:  F1LE2 

Variable  Value  Meaning 


! 


QL2L(K) 


Initial  lateral  inflow  discharge  at  all  reaches  containing  lateral  inflow 


Statement  19 

Variables:  QCHECI(I) 
Format:  FIO.O 
File  read:  FILE2 
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Variable 


Value  Meaning 


QCHECI(I)  Initial  discharge  of  dam  that  has  a  rating  curve  as  input 


There  is  one  card  no.  19  for  every  dam  that  has  a  rating  curve  as  input. 


Statement  20 

Variables:  NLEVEE,  ILEVEE(I),  1=1,  NLEVEE 
Format:  1615 
File  read:  FILE2 

Variable  Value  Meaning 


NLEVEE  0  Number  of  reaches  with  levees 

ILEVEE(I)  Upstream  nodal  points  of  reaches  with  levees 


Statement  21 

Variables:  ELEVEE(K),  K=l,  NLEVEE 
Format:  8F10.0 
File  read:  FILE4 

Variable  Value  Meaning 


ELEVEE(K)  Average  elevation  of  the  top  of  the  levee  along  this  reach  in  feet 


If  NLEVEE  =  0,  this  card  is  skipped. 


Statement  22 

Variables:  RANGE(I),  XL(I),  DUMMY.  ZF(I),  Z(I),  BETA(I) 
Format:  A8,  2X,  7F10.0 
File  read:  FILES 


Variable  Value  Meaning 


RANGE(I) 

Description  of  the  /*'’  nodal  point 

XL(I) 

River  mileage  of  nodal  point 

DUMMY 

Space  for  top  width;  can  be  left  blank 

ZF(I) 

Top  bank  elevation  for  nodal  point  in  feet 

Z(I) 

Bed  elevation  of  nodal  point  in  feet 

BETA(I) 

Momentum  correction  factor 

Statement  23 

Variables:  H1(U),  AI(I,J),  TI(U),  RNI(LJ) 
Format:  4F10.0 
File  read:  FILES 
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Variable 


Value  Meaning 


AI(I,J) 

TI(I,J) 

RNI(I,J) 


elevation  of  channel  geometry  tables  for  point  I 
Flow  area  corresponding  to  elevation 
Top  width  corresponding  to  elevation 
Manning’s  n  corresponding  to  elevation 


Card  no.  23  is  repeated  ISTAGE  times. 


Statement  24 

Variables:  HF{I,J),  AHd.J),  DUMMY,  RNIFP(1.J) 
Format:  4F10.0 
File  read:  FILES 


Variable  Value 

Meaning 

HF(I,J) 

AFI(I,J) 

DUMMY 

RNIFP(U) 

elevation  of  the  channel  flood  plain  geometry  table 

Row  area  corresponding  to  elevation  HF(1  J) 

Top  width  corresponding  to  elevation  HF(1  J),  can  leave  blank 
Manning’s  n  corresponding  to  elevation  HF(I  J) 

Card  no.  24  is  repeated  IFPLN  times. 

Cards  no.  22,  23  and  24  are  repeated  for  every  nodal  point. 

Statement  25 

Variables:  NSEG 

Format:  IS 

File  read:  FILE4 

Varirll^  Value 

Meaning 

NSEG  6 

Number  of  linear  segments  approximating  the  rating  curve 

Statement  26 

Variables:  QRC(I),  HRC(I) 
Format:  2F10.0 

File  read:  nLE4 

Variable  Value 

Meaning 

QRC(rr 

HRC(l) 

Discharge  at  the  end  of  the  /’*’  linear  segment 

Elevation  of  water  surface  corresponding  to  the  end  of  the  segment 

Card  no.  26  is  repeated  NSEG  times. 
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Statement  27 

Variables:  KRC(I),  QLIMIT(I),  QCHECKO(I),QDRCF(I),  QDRC(K,  I),  HDRC(K,I),  K=l,  KRC(I) 
Format;  15,  7F10.0/8F10.0 
File  read;  FILE4 

Variable  Value  Meaning 


KRC(I) 

QLIMIT(I) 

QCHECKO(I) 

QDRCF(I) 

QDRC(K,I) 

HDRC(K,I) 


Number  of  entries  in  rating  curve  table  at  the  /**’  dam 

Discharge  below  which  a  fixed  water  surface  elevation  is  prescribed 

Initial  discharge  of  the  /'*'  dam 

Discharge  above  which  the  falling  portion  of  the  rating  curve  will  he  used 

if  discharge  is  decreasing 

/f**’  discharge  in  the  rating  curve  table  of  dam  / 

Water  surface  elevation  corresponding  to  QDRC(K,I) 


Card  no.  27  is  repeated  for  each  dam  with  NVARY(I)  =  2. 


Statement  28 

Variables:  J,  QL2(K),  K=l,  NRCH 
Format:  15,  7F10.0/8F10.0 
File  read;  FILEl 

Variable  Value  Meaning 


J  24  Number  of  time  steps  before  new  lateral  inflows  will  be  input 

QL2(K)  Lateral  inflow  in  cfs 


Statement  29 

Variables:  S,  J 
Format;  FI 0.0, 15 
File  read:  FILEl 

Variable  Value  Meaning 


S  New  downstream  boundary  condition 

J  Number  of  time  steps  before  a  new  boundary  condition  will  be  input 


Statement  30  '< 

Variables:  GRAVIT,  ROW,  ROI,  SPHT,  XLATEN 
Format:  5F10.0 
File  read:  FILE6 

Variable  Value  Meaning 


GRAVIT  9.8 1  Acceleration  of  gravity  in  m/s^ 

ROW  1000  Density  of  water  in  kg/m^ 
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ROI  916  Density  of  ice  in  kg/tn-^ 

SPHT  4185.5  Specific  heat  of  water 

XLATEN  334000  Latent  heat  of  melting  of  ice 


Statement  31 

Variables:  ALPCOV,  PORFRA,  PORINI,  VEROS.  VDEPOS.  VBUOY 
Format:  6F10.0 
File  read:  FILE6 


Variable 

Value 

Meaning 

ALPCOV 

0.85 

Fraction  of  concentration  going  for  ice-cover  formation 

PORFRA 

0.6 

Porosity  of  frazil  ice 

PORINI 

0.2 

Porosity  of  initial  ice  cover 

VEROS 

0.6 

Velocity  of  erosion  in  m/s 

VDEPOS 

0.5 

Velocity  of  deposition  in  m/s 

VBUOY 

0.001 

Buoyant  velocity  in  m/s 

Statement  32 

Variables:  XKS,  XKI,  ES,  TOLNUL 

Format:  4F10.0 

File  read:  FILE6 

Variable 

Value 

Meaning 

XKS 

0.25 

Thermal  conductivity  of  snow 

XKI 

2.24 

Thermal  conductivity  of  ice 

ES 

0.6 

Porosity  of  snow 

TOLNUL  0.00001 

Small  number  used  strictly  for  programming  purposes 

Statement  33 

Variables:  XMU,  XK2,  COHE,  THIBLK,  HWA,  SIGMA 

Format:  6F10.0 

File  read:  FILE6 

Variable 

Value 

Meaning 

XMU 

1.28 

p  used  in  Pariset  and  Hausser’s  wide  jam  equation 

XK2 

3.0 

Coefficient  of  passive  stress  for  granular  ice 

COHE 

100.0 

Bank  cohesion  used  in  the  wide  jam  equations 

THIBLK 

0.08 

Initial  thickness  of  ice  block  in  m 

HWA 

23.9 

Water  to  air  heat  transfer  coefficient 

SIGMA 

81550 

Compressive  strength  of  solid  ice  crust  formed  due  to  cooling  from  top 
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Statement  34 

Variables:  LFnM(lRCH(I)),  1=1,  NR 
Format:  1615 
File  read:  FILE6 

Variable  Value  Meaning 


LFT1M(IRCH(I))  1  Lateral  inflow  will  have  temperature  assigned 

0  Lateral  inflow  will  not  have  temperature  assigned 


Statement  35 

Variables:  KDUM,  IBOOM(I).  DXICM(I).  RICEM(I),  FRPM(I).  FRCM{I) 
Format:  215, 4F10.0 
File  read:  FILE6 

Variable  Value  Meaning 


KDUM 

Number  of  the  reach 

IBOOM(I) 

1 

Ice  bridging  possible  at  the  end  of  the  reach 

0 

Ice  bridging  not  possible  at  the  end  of  the  reach 

DXICM(I) 

Initial  length  of  the  ice  cover  in  reach  / 

RICEM(I) 

Ice  cover  roughness  on  reach  / 

FRPM(I) 

0.06 

Critical  Froude  number  for  undenuming 

FRCM(I) 

0.09 

Critical  Froude  number  for  progression 

Card  no.  35  is  repeated  for  every  reach. 


Statement  36 

Variables:  NXT 
Format:  15 
File  read:  FILE6 

Variable  Value  Meaning 


NXT  13  No.  of  nodal  points  where  air  temperature  prescribed  as  boundary 

condition 


Statement  37 

Variables:  NTA(I),  1=1,  NXT) 
Format:  1615 
File  read:  FILE6 

Variable _ Value  Meaning 


NTA(I) 


Number  of  nodal  point  where  air  temperature  is  prescribed 
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Statement  38 

Variables;  NTW(I),  1= 
Format:  1615 

File  read:  FILE6 

LNXTW 

Variable  Value 

Meaning 

NTW(I) 

Number  of  nodal  point  where  wa'er  temperature  is  prescribed 

Statement  39 

Variables:  IDUMM 
Format;  15 

File  read:  FILES 

Variable  Value 

Meaning 

IDUMM  0 

1 

No  initial  ice  cover 

Existing  initial  ice  cover 

Blank  card  follows. 

Statement  40 

Variables:  IDUM,  DXICM(I),  TICEM(I),  TFRAM(I).  TSUBM(I),  THINIM(I) 

Format:  15,  6F12.5 

File  read:  FILES 

Variable  Value 

Meaning 

IDUM 

DXICM(I) 

TICEM(I) 

TFRAM(I) 

TSUBM(I) 

THINIM(I) 

Reach  number 

Ice  cover  length  coefficient  on  reach  / 

Solid  ice  cover  thickness  in  reach  /  in  m 

Frazil  ice  thickness  in  reach  /  in  m 

Submerged  ice  cover  thickness  in  reach  /  in  m 

Initial  ice  cover  thickness  in  reach  /  in  m 

There  is  one  card  no.  S9  for  every  reach. 

If  IDUMM  =  0  (card  no.  S8),  there  will  be  no  card  no.  S9. 

A  blank  card  follows. 

Statement  41 

Variables:  TWS(I),  1= 
Format:  10F7.0 

File  read:  FILES 

1,NX 

Variable  Value 

Meaning 

TWS(I) 

Initial  water  temperature  at  every  nodal  point 

Two  blank  cards  follow. 
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Statement  42 

Variables:  TAS(I),  I=:1,NX 
Format:  10F7.0 
File  read:  FILES 


Variable  Value 

Meaning 

TAS(I) 

Initial  air  temperature  at  every  nodal  point 

There  is  one  blank  card  before  card  no.  4 1 . 

Statement  43 

Variables:  TWXA(I),  1= 
Format:  8F10.0 

File  repd;  HI  E6 

d,NJUNC+l 

Variable  Value 

Meaning 

TWXA(I) 

Boundary  water  temperature  at  the  upstream  end  of  /**’  tributary 

Statement  44 

Variables:  TLIM 

Format:  FI 0.0 

File  read:  FILE6 

Variable  Value 

Meaning 

TLIM(I) 

Temperature  of  the  lateral  inflow 

Card  no.  43.  is  repeated  for  every  lateral  inflow  that  had  LFTlM(l)  ^  0. 

Statement  45 

Variables:  TAS(NTA(I)),  1=1,  NXT 

Format:  8F10.0 

File  read:  FILE6 

Variable  Value 

Meaning 

TAS(NTA(I)) 

New  boundary  air  temperature  at  point  NTA(I) 

Cards  no.  44,  42  and  43  are  repeated  for  every  day  of  computation. 
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Output 

Only  printed  output  can  be  obtained  as  a  result  of  the  program.  Several  plotting  programs  are 
available  to  represent  these  results  in  graphic  mode  when  needed.  RESl  .DAT  contains  all  important 
parameters  and  descriptions  of  the  system,  followed  by  the  time  history  of  the  water  surface  elevation, 
flow  discharge  and  velocity,  and  bed  elevation  at  each  nodal  point  requested.  RES2.DAT  contains  the 
time  history  of  the  water  surface  elevation  and  flow  discharges  for  all  nodal  points,  the  lateral  inflows 
and  the  discharges  for  dams  with  rating  curves  as  input.  This  file  can  be  used  as  input  FILE2.DAT 
without  any  changes.  RES3.DAT  indicates  whether  the  run  was  successful  or  not.  RES4.DAT  contains 
water  and  air  temperatures  for  every  nodal  point  along  the  river,  followed  by  ice  cover  conditions. 
These  include  the  reach  number,  the  fraction  of  reach  under  an  ice  cover,  the  thickness  of  the  solid  ice 
cover,  the  thickness  of  frazil  deposition,  and  the  submerged  and  initial  ice  cover  thickness.  All 
thickness  data  are  given  in  meters.  If  an  ice  cover  is  not  present  on  the  river,  information  about  ice  cover 
conditions  is  omitted.  This  file  can  be  used  as  input  FILE3.DAT. 

Files  RES1.DAT,  RES2.DAT  and  RES3.DAT  are  exactly  the  same  as  the  output  files  from  the 
original  hydraulics  model,  so  they  will  not  be  discussed  further.  Detailed  descriptions  of  these  files  are 
provided  by  Johnson  ( 1 982).  A  sample  output  of  file  RES4.DAT  is  given  on  the  following  pages.  Since 
the  results  of  the  computations  are  given  for  each  reach  separately,  it  is  not  easy  to  interpret  them  for 
every  pool.  A  small  FORTRAN  program  is  provided  to  rearrange  these  data  in  a  simpler  pattern.  The 
listing  of  this  program,  together  with  its  output,  is  given  after  the  sample  of  output  file  RES4.DAT. 

Limitations 

RICEOH  is  a  one-dimensional  unsteady-flow  thermal  and  ice  model  that  can  be  applied  to  a  general 
system  of  river  channels  containing  locks  and  dams.  However,  there  are  some  limitations  to  its 
applicability.  This  model  can  be  applied  only  to  a  simply  connected  system,  i.e.  closed  loops  within 
the  system  cannot  be  handled.  An  additional  limitation  on  the  physical  system  is  that  there  can  be  only 
one  downstream  boundary. 

In  its  present  form,  there  is  some  limitation  on  the  specification  of  boundary  conditions.  At  an 
upstream  boundary,  only  flow  discharges  can  be  prescribed;  at  a  downstream  boundary,  either  a  rating 
curve  or  water  surface  elevations  may  be  specified.  In  general,  one  can  specify  elevations  at  the 
upstream  boundary  and  discharges  at  the  downstream  boundary,  but  some  additional  modifications 
on  the  current  version  will  be  required. 

The  major  limitation  on  the  thermal  and  ice  routines  is  that  negative  velocities  cannot  be  handled. 
If  any  of  the  branches  at  any  time  flow  upstream,  the  temperature  computation  cannot  be  performed. 
Since  the  time  period  of  this  simulation  is  mainly  in  the  winter  (the  low-flow  season),  this  limitation 
is  not  important. 

Variable  names 

In  this  section  definitions  of  input  as  well  as  other  variables  in  the  program  are  tabulated  in  Tables 
A 1  and  A2. 
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Table  Al.  Definitions  and  typical  values  of  input  variables. 


Variable 

Value 

Definition 

Cross-sectional  area  at  elevation  Hl(ij)  (ft');  i  =  1,2,..  .NY: ;  -  1 
2,...ISTAGE. 

AFI(iJ) 

Cross-sectional  area  of  the  flood  plain  at  elevation;  HF(iJ)  (ft');  / 
=  \  ,2,...NXJ=  \,2,.../FPLN. 

AL(i) 

1.0 

Velocity  head  correction  factor  associated  with  junction  of  up¬ 
stream  main  river  and  downstream  main  river;  /  -  ’ ,  2,  NJUNC 

*  21. 

ALU) 

1.0 

Velocity  head  correction  factor  associated  with  junction  of  tribu¬ 
tary  and  downstream  main  river;  /  =  2,  4,. .  .NJUNC  *  2-1 . 

ALPCOV 

0.90 

Fraction  of  total  ice  flow  going  into  cover  formation  (0-1). 

BETA(i) 

1.0 

Correction  factor  in  momentum  equation; ;  =  1,  2, ...NX. 

COHE 

100.0 

Bank  cohesion  used  in  the  wide  jam  equations  (kg  m“-). 

DXICMU) 

Fraction  of  the  length  of  a  river  reach  covered  by  ice  (0.0-1 .0); 

;■  =  1,  2,... NR.  Before  the  formation  of  ice  cover,  as  input,  set  this 
vanable  to  TOLNUL  at  all  points  where  ice  cover  progression  is 
possible. 

DUM 

Dummy  variable. 

ELEVEEU) 

Average  elevation  of  the  top  of  the  levee  along  reach  /;  /  =  1, 
2,...NLEVEEifX). 

ES 

0.6 

Porosity  of  snow  (0.0- 1.0). 

FRCMU) 

Critical  Froude  number  above  which  ice  cover  progression  is  not 
possible  in  reach  /;  i=  1,2,..  .NX. 

FRPMii) 

Froude  number  below  which  particle  juxtraposition  takes  place  in 
the  river  reach  /;  /  =  1,2, ...NX. 

G 

32.1614 

Gravitational  acceleration  (ft  s“^). 

GRAVIT 

9.8 

Gravitational  acceleration  (m  s”-). 

Initial  water  surface  elevation  at  each  net  point  (ft);  i=  1,2,..  .NX. 

HDRC(j,i) 

Water-surface  elevation  corresponding  to  QDRC(j,i)  (ft);  7  =1, 
2,...KRC(,i). 

mij) 

Elevation  of  the  flood  plain  geometry  (ft);  /  =  1,  2,...NX\  7  =1, 
2,...IFPLN. 

mij) 

Elevation  of  channel  geometry  table  i  =  1,  2, ...NX-,  7=1, 
2,..  1ST  AGE. 

HRC{i) 

Elevation  of  water  surface  corresponding  to  the  end  of  the  linear 
segment  i  (ft);  i  =  1,2...  .NSEG. 

HSETO(i) 

Elevation  maintained  by  dam  (ft). 

HWA 

23.9 

Water-to-air  heat  transfer  coefficient  (W  m““  °C'’). 

IBACK 

1 

Normal  output  (if  IBACK  =  0  less  printed  output). 

IBOOM(i) 

Flag  to  indicate  the  presence  of  an  ice  boom,  dam  or  ice  control 
structure  at  the  end  of  reach  i.  If  =  1 ,  ice  cover  can  progress  starting 
from  the  downstream  end  of  reach  v,  i=  1,2...  .NX. 

IBC(i) 

-1 

Rating  curve  will  be  used  if  this  branch  contains  an  outer  down¬ 
stream  boundary. 

0 

This  is  an  interior  branch,  or  discharge  is  input  at  its  upstream 
boundary. 

1 

Elevations  will  be  input  if  this  branch  contains  an  outer  downstream 
boundary. 
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Table  A1  (cont’d). 


Variable 

Value 

Definition 

IBRNCH{i,l) 

First  net  point  on  branch. 

IBRNCHU,2) 

Last  net  point  on  branch. 

IBRS(i) 

0 

Branch  is  an  interior  branch. 

1 

Branch  has  an  upstream  outer  boundary. 

ID 

Dummy  argument. 

lEDYHD 

0 

Eddy  head  loss  coefficients  set  to  zero  (if  lEDYHD  =  1 ,  eddy  head 
loss  coefficients  are  read  in). 

IFPLN 

3 

Number  of  entries  in  flood  plain  tables. 

IGEOM 

1 

Geometry  tables  will  not  be  printed  (if  IGEOM  =  0,  they  will  be 
printed). 

IHSET 

24 

Number  of  time  steps  before  a  new  water  surface  elevation  to  be 
maintained  by  dam  is  read. 

IJUNCU,\) 

Number  of  upstream  branch  on  main  river  before  junction 
/■=  2,... NJUNC. 

IJUNC(i,2) 

Number  of  tributary  branch  i  =  1,2,..  .NJUNC. 

IJUNCU,3) 

Number  of  downstream  branch  on  main  river  after  junction 
/=  1,2,... NJUNC. 

ILEVEEii) 

Upstream  net  points  of  reaches  with  levees. 

ILUG 

90 

Unit  from  which  geometry  tables  will  be  read. 

INTVD 

24 

Print  interval  for  particular  days. 

INTVG 

24 

Major  print  interval. 

INTVP 

24 

Interval  for  placing  points  in  plot  file. 

IPLT 

0 

No  plots. 

IPRINT 

0 

Limited  output  (if  IPRINT  =  1 ,  detailed  output). 

IQCK 

24 

Number  of  time  steps  before  a  new  discharge  will  be  input. 

IRCH(i) 

Numbers  of  the  reaches  containing  lateral  inflow  i  =  1 , 2,. .  .NRCH. 

ISTAGE 

5 

Number  of  entries  in  channel  geometry  tables. 

J 

24 

Number  of  time  steps  before  new  lateral  inflows  will  be  input. 

KRC(i) 

15 

Number  of  entries  in  rating  curve  table  at  the  dam  i. 

LFTIM(i) 

Indicator  of  thermal  importance  of  lateral  inflow;  if  LFTIM{i)  =  1, 
lateral  inflow  i  is  thermal  inflow. 

NBRS 

9 

Total  number  of  branches. 

NC 

185 

Total  number  of  net  points  minus  one. 

NDAMS 

13 

Total  number  of  dams. 

NJUNC 

4 

Total  number  of  junctions. 

NLii) 

Net  point  immediately  upstream  of  dam. 

NLEVEE 

0 

Upstream  net  points  of  reaches  with  levees. 

NOXS 

0 

Number  of  stations  at  which  plots  are  desired. 

NPRINT(i) 

Net  point  numbers  at  which  output  is  desired;  i  =  1,  2, ...NSTAT. 

NRCH 

38 

Total  number  of  reaches  containing  lateral  inflows. 

NSEG 

6 

Number  of  linear  segments  approximating  the  rating  curve. 

NSTAT 

37 

Number  of  net  points  at  which  output  is  desired. 

NTA(i) 

Net  point  where  boundary  air  temperature  is  prescribed. 

NTILF 

1 

Number  of  lateral  inflows  that  are  thermal  inflows. 

NVARY(i) 

0 

Normal  dam. 
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Table  A1  (cont’d).  DeFinitions  and  typical  values  of  input  variables. 


Variable 

Value 

Definition 

1 

Time-varying  elevations  of  pool  will  be  input. 

2 

Rating  curve  will  be  input  for  this  dam. 

NXMAIN 

134 

Last  net  point  on  main  river. 

NXT 

13 

Total  number  of  points  where  boundary  air  temperature  is  pre¬ 
scribed. 

PORFRA 

0.4 

Porosity  of  frazil  ice  accumulation  (0.0-1. 0). 

PORINI 

0.2 

Porosity  of  initial  ice  cover  (0.0- 1.0). 

G(/,l) 

Initial  discharge  at  each  net  point  (ft^  s“’); ;  =  1, 2, ...NX. 

Q(i,2) 

New  boundary  discharge  at  upstream  end  of  tributaries  and  main 
stream  (ft^  s“*). 

QCHECO(i) 

Initial  discharge  of  dam  /'  (ft^  s“'). 

QDRCFii) 

Discharge  above  which  the  falling  portion  of  the  rating  curve  for 
dam  i  will  be  used  (ft^  s"'). 

QDRCijJ) 

Discharge  in  the  rating  curve  table  for  dam  i  (ft^  s“');  7=1, 
2,...KRC{i). 

QL2(i) 

Lateral  inflow  (ft^  s“');  /  =  1,2,.. .NRCH. 

QLIMITU) 

Discharge  below  which  a  fixed  water-surface  elevation  for  dam  i  is 
prescribed  (ft^  s~*). 

QRcm 

Discharge  at  the  end  of  the  linear  segment  i  (fp  s”');  /  =  1, 
2,...NSEG. 

RANGEii) 

Description  of  net  point  i;  /  =  1,  2,... NX. 

RICEMU) 

Ice  roughness  for  reach  i;  i  =  1,2,..  .NX. 

RNI(i) 

Manning’s  roughness  coefficient  at  the  elevation  HI{ij)-, 

1  =  1 , 2,. .  .NX;  7  =  1 , 2 . 1ST  AGE. 

RNlFPii) 

Manning’s  roughness  coefficient  at  the  elevation  HF{iJ), 

/=  1,2 . NX,j=  \,2,...IFPLN. 

ROI 

916.0 

Density  of  ice  (kg  m"^). 

ROW 

1000.0 

Density  of  water  (kg  m"^). 

SIGMA 

0.408  X  105 

Compressive  strength  of  solid  ice  crust  formed  due  to  cooling  from 
top  (kg  m“^). 

SPHT 

4.1855  X  10^ 

Specific  heat  of  water  (J  kg"’  °C“’). 

TAS(i) 

Air  temperature  at  net  point  i  (°C);  (  =  1,  2, ...NX. 

TDAM(i) 

Description  of  the  dam. 

TFRAM(i) 

Thickness  of  frazil  deposition  on  reach  r  (m);  i  =  1, 2,... NX;  will  be 
input  only  if  IDUMM  =  1 . 

THIBLK 

0.08 

Thickness  of  a  floating  ice  floe  (m). 

THINIMii) 

Thickness  of  initial  ice  cover  on  reach  i  (m);  i  =  1, 2,...A/X;  will  be 
input  only  if  IDUMM  =  1 . 

TIUJ) 

Top  width  in  geometry  table  /  =  1, 2, ...NX; 7  =1,  2,...ISTAGE. 

TITLE 

Description  of  run. 

TICEM(i) 

Thickness  of  solid  ice  cover  on  reach  /;  /  =  1,2,..  .NX;  will  be  input 
only  if  IDUMM  =  1  (m). 

TLIM(i) 

Water  temperature  of  lateral  inflow  that  has  LFTIMfi)  =  1  (°C). 

TOLNUL 

0.00001 

Very  small  positive  non-zero  number  used  in  programming. 

TOTALT 

Number  of  days  of  computations. 
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Table  A1  (cont’d). 


Variable 

Value 

Definition 

TSTEP 

TSUBMii) 

3600 

Time  step  (s). 

Thickness  of  submerged  ice  cover  on  reach  /  (m);  /  =  1 , 2,. .  .NX',  will 
be  input  only  if  IDUMM  =  1 . 

TT 

TWS{i) 

86400.0 

Length  of  time  step  used  in  the  simulation  (s). 

Initial  water  temperature  at  net  point  i  (°C);  i  =  1,  2, ...NX. 

VBUOY 

0.001 

Buoyancy  velocity  of  frazil  ice  particles  used  in  determining  the 
amount  of  deposition  (m  s”*). 

VDEPOS 

0.6 

Velocity  of  deposition,  which  is  the  velocity  below  which  deposi¬ 
tion  of  frazil  ice  under  the  ice  cover  is  possible  (m  s“'). 

VEROS 

0.7 

Velocity  of  erosion,  which  is  the  velocity  above  which  erosion  of 
frazil  ice  under  the  ice  cover  is  possible  (m  s“*). 

XKI 

2.24 

Thermal  conductivity  of  ice  (W  m"’  “C'"'). 

XKS 

0.25 

Thermal  conductivity  of  snow  (W  m“'  °C“'). 

XK2 

XL{i) 

3.00 

Coefficient  of  passive  stress  for  granular  ice;  XK2  =  tan-  (n/4  -1-  (j)/2). 
River  mileage  of  net  point  1;  /  =  1 , 2,. .  .NX.',  tributary  mileage  is  zero 
at  junction  (miles). 

XLATEN 

3.34  X  10^ 

Latent  heat  of  melting  of  ice  (J  kg“'). 

XMU 

ZFU) 

Z(i) 

1.28 

|i  used  in  Pariset  and  Hausser’s  wide  Jam  equation. 

Top  bank  elevation  for  net  point  1  (ft);  /  =  1,2,..  .NR. 

Bed  elevation  for  net  point  i  (ft);  i  =  1,  2, ...NR. 

Table  A2.  Derinitions  of  \  ariables  other  than  input  variables  used  in  the  program. 


Variable 


Definition 


Aii,\) 

A(/,2) 

ABU) 

AD[i) 

AF{i) 

AU(i) 

AXyCl.l) 

AXy(/,2) 

BMAN(i) 

CO(iJ) 

CSI 

Eli) 

FRCONV 


Flow  area  as  calculated  in  hydraulic  part  during  the  current  time  step  (ft“);  i  = 
\,2,...NX. 

Flow  area  as  calculated  in  hydraulic  part  during  the  previous  time  step  (ft“);  i 
=  \,2,...NX. 

Average  effective  flow  area  of  reach  i.  This  area  is  used  in  water  temperature 
and  ice  computations  (m^);  t  =  1, 2,... NR. 

Downstream  effective  flow  area  of  reach  i  =  1 ,  2,. .  .NR  (m^). 

Flood  plain  area  as  calculated  in  hydraulic  part  (ft“);  i  =  1, 2,... NX. 

Upstream  effective  flow  area  of  reach  i;  /  =  1,  2, ...NR  (m‘). 

Derivative  of  flow  area  along  the  channel  when  flow  depth  is  kept  constant 
during  the  current  time  step;  /  =  1,2,..  .NX. 

Derivative  of  flow  area  along  the  channel  when  flow  depth  is  kept  constant 
during  the  previous  time  step;  i  =  \,2,...NX. 

Manning’s  roughness  coefficient  for  reach  i;  i  =  1,  2, ...NR. 

Coefficients  of  the  two  governing  equations;  i=  1,2,..  .NX*2;  j  =  1,2,3,  4. 
Parts  oiCO{ij)\l=  1,2....  13. 

Right  side  coefficients  of  the  two  governing  equations;  i  =  1,  2,...NX*2. 
Factor  to  convert  temperature  to  ice  concentration. 
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Table  A2.  Definitions  of  variables  other  than  input  variables  used  in  the  program. 


Variable 


Definition 


lOAREA(iJ) 


D(i,l) 

D(/.2) 

DHH 

DHT 

DXIC(i) 

HHD{i) 

HHU{i) 

IBOO(i) 

/PROG(i) 

JEP 

JSP 

NB 

NOAREA 

NR 

QBii) 

QD(i) 

QLAT(i) 

QLATKi) 

QUO) 

R/CE(i) 

RMUNO) 

RMUNKi) 

S(0 

SF(i,\) 

SF0,2) 

m) 

ni,2) 

TB(i) 

TOO) 

TFRAO) 

THINHi) 


Variable  to  designate  the  extent  of  open-water  patches  in  an  ice-covered  riven 
j  =  counter  for  the  open  water  are  starting  from  upstream  (<  NR)-,  i  =  1  gives  the 
starting  node  number  for open  water  area; ;  =  2  gives  the  ending  node  number 
of /^open  water  area; « =  3  gives  the  last  node  up  to  which  a  particle  in  this  open- 
water  area  will  travel  during  the  time  step. 

Depth  of  the  flow  in  net  point  /  during  the  current  time  step  (ft);  i=  1,2,..  .NX. 
Depth  of  the  flow  in  net  point  /  during  the  previous  time  step  (ft); ;  =  1 , 2,. .  .NX. 
Value  of  5//  used  when  calculating  deri  vati  ves  of  geometric  values  with  respect 
to  the  stage  (ft). 

Twice  the  value  of  DHH,  used  in  calculating  derivatives  of  geometric  values 
with  respect  to  the  stage  (ft). 

Fraction  of  the  length  of  a  river  reach  covered  by  ice;  i=  1,2,..  .NR. 
Downstream  water  level  of  reach  /  (m);  /  =  1,2,..  .NR. 

Upstream  water  level  of  reach  /  (m);  /  =  1,2,..  .NR. 

Flag  to  indicate  the  presence  of  an  ice  boom,  dam  or  an  ice  control  structure  at 
the  end  of  the  reach  /;  i  =  \  ,2,...NR. 

Flag  to  indicate  whether  progression  is  possible  at  reach  /  after  checking  jam 
conditions;  0  =  no;  I  =  yes;  i  =  1.2,..  .NR. 

Variable  used  in  controlling  the  time  step. 

Variable  used  in  controlling  the  time  step. 

Number  of  nodal  points  in  branch  while  performing  thermal  computations. 
Number  of  open-water  patches  in  the  entire  river;  see  lOAREA. 

Number  of  reaches  in  branch  while  performing  thermal  computations. 
Average  flow  rate  of  reach  /  (m^  s"');  <  =  1,2,..  .NR. 

Discharge  at  downstream  end  of  reach  i  (m^  s"’);  /  =  1,  2, ...NR. 

Lateral  inflow  discharge  of  reach  i  (m^s“’);  /  =  1,  2,.  ..NX. 

Lateral  inflow  discharge  of  reach  i  (m^  s"*);  /  =  1,  2, ...NR. 

Discharge  at  upstream  end  of  reach  i  (m^  s"');  i  =  1,  2,. ..NR. 

Ice  roughness  (Manning)  of  a  reach ;;  /  =  1  ,  2, ...NR. 

Manning’s  roughness  coefficient  of  a  point  /;  /  =  1,  2, ...NX. 

Fraction  of  an  ice  cover  during  the  previous  time  step  at  reach  /. 

Distance  from  node  I  to  termination  point  of  t Lagrangian  moving  point  at  the 
end  of  the  time  step  (m);  /  =  1,  2,...NB. 

Frictional  slope  at  net  point  i  during  the  current  time  step;  /  =  1,2,..  .NX. 
Frictional  slope  at  net  point  i  during  the  previous  time  step;  t  =  1,  2, ...NX. 
Channel  top  width  at  net  point  /  during  the  current  time  step  (ft);  i=  1,2,..  .NX. 
Channel  top  width  at  net  point  i during  the  previous  time  step  (ft);  t  =  1 , 2,. .  .NX. 
Average  top  width  of  river  section  i  used  in  the  ice  computations  (m);  /  =  1, 

2.. .. NR. 

Downstream  top  width  of  river  section  i  used  in  the  ice  computations  (m);  i  = 

1.2.. .. NR. 

Thickness  of  frazil  ice  in  reach  i  (m);  i  =  1,  2, ...NR. 

Thickness  of  initial  ice  cover  in  reach  i  (m);  i  =  1,  2. ...NR. 
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Table  A2  (cont’d). 


Variable 


Definition 


THPROii) 


TlCE{i) 

TLIU) 

TSUBii) 

TUU) 

TW(i) 


TWXAU) 

TWOm) 


U(i,]) 

U(i,2) 

UBii) 

WP{i) 

WPIU) 

X(i) 

ZD{i) 

ZU(i) 


Thickness  of  initial  ice  cover  of  reach  i  computed  using  jam  equations.  This 
thickness  is  used  to  compute  the  length  of  progression  for  a  given  volume  of  ice 
(m);  /  =  1,  2,...NR. 

Sum  of  black  and  white  ice  thicknesses  in  reach  /  (m);  /  =  1,2,..  .NR. 
Temperature  of  the  lateral  inflow  /  during  the  thermal  computations. 
Submerged  thicknesses  of  the  ice  cover  in  reach  /  (m);  /  =  1 ,  2, . .  .NR. 
Upstream  top  width  of  river  section  i  used  in  the  ice  computations  (m);  /'  =  1 , 

2.. .. NR. 

Cross-sectional  average  water  temperature  at  node  i  (°C).  The  same  variable  is 
used  to  express  cross-sectional  average  frazil  ice  concentration  at  node  i  as  a 
temperature.  When  71V’(/)  is  negative,  concentration  is  obtained  by  multiplying 
it  by  FRCONV. 

Water  temperature  of  node  1  at  next  time  step  (°C). 

Individual  frazil  concentration  distributions  of  each  of  the  open-water  areas; 

1  =  node  number;  j  =  open-water  patch  counter;  /  =  1,  2,...NB;  7  =1, 

2. . ..NOAREA.  This  variable  is  used  in  subroutine  OPENLAGR. 

Flow  velocity  at  net  point  i  during  the  current  time  step  (ft^s*’);  /  =  1 , 2, ...NX. 
Flow  velocity  at  net  point  /  during  the  previous  time  step  (ft^s"' );  /  =  1 , 2,. .  .NX. 
Average  river  flow  velocity  in  reach  i  (m  s"’);  i  =  1, 2, ...NR. 

Wetted  perimeter  of  node  /  (ft);  i  =  1,2,..  .NX. 

Wetted  perimeter  of  node  i  during  the  previous  time  step  (ft);  i  =  1,2,..  .NX. 
Distance  from  node  1  to  node  i  (m);  /  =  1  ,  2,. .  .NR. 

Height  of  the  reference  elevation  at  the  downstream  end  of  the  reach  i  (m);  i  = 
I,  2,... NR. 

Height  of  the  reference  elevation  at  the  upstream  end  of  the  reach  i  (m);  /  =  1 , 

2.. .. NR. 
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